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ABSTRACT 


The performance of Loran-C receivers is severely limited by 
contamination introduced by multipath propagation effects. Echo removal 
by discrete generalized linear filtering is considered as a method of 
reducing the effects of the contamination. The signals can be described 
as consisting of repeated, known, composite signals with short echo 
arrival times, in noise. Since the signals are repeated, and stable, 
they can be averaged to increase the effective signal to noise ratio. 
The homomorphic deconvolution technique is applied to simulated signals 
and a determination of the averaging time needed for satisfactory 
performance is made, 

Aside from the Loran-C application, the report presents a detailed 
Study of the computational aspects involved in deriving the complex 
cepstrum, Since the signal is known, the effects of noise and aliasing 
can be precisely identified. In particular, conventional phase 
unwrapping methods are questioned for certain classes of signals. 
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Chapter 1 


INTRODUCTION 


The digital signal processing technique referred to a Homomorphic 
Deconvolution has met with some success, in a number a diverse appli- 
cations, in suppressing signal contamination introduced by a multi- 
path propagation medium. Another possible application is to Loran-C, 
a low frequency electronic navigation system, in which the performance 
is severely limited by the effects of such contamination. It is the 
purpose of this thesis to study the applicability of the technique to 
various aspects of Loran-C signal processing, 

In this Chapter, a brief description of the Loran-C system and 
Homomorphic Deconvolution is presented. In Chapter 2, the task of the 
homomorphic filter in the Loran-C case is formulated. In Chapter 3, 
the filtering technique is developed and the results of simulated 
filtering in the noiseless case are presented. Chapter 4 presents the 
theory and a Simulated example of the effects of signal distortion, 
Chapter 5 includes a theory of the effect of white gaussian noise on 
the processing and includes examples to verify the theory. In 
Chapter 6, the effects of the various assumptions and simplifications 
made are discussed and the results of the preceeding chapters are 
examined to determine the applicability of the method, The examples 
presented are simulations performed on an IBM 370 general purpose 


computer. 





1,1 Loran-C System Description b 

Loran-C is an electronic navigation system which belongs to a 
broad class of systems referred to as hyperbolic, a characterization 
based on straightforward geometric considerations. If radio signals 
are simultaneously transmitted from two different locations, they will, 
in general, arrive at some arbitrary third location at two different 
times. To the extent that the speeds of propagation along the two 
paths are equal and known, a measurement of the difference in times 
of arrival provides a corresponding measure of the difference in 
distances involved. The locus of points with a constant difference 
in distance from two reference points is a hyperbola. A receiver 
which determines this time difference and identifies the proper 
reference points provides information sufficient to define a navi- 
gational line of position. Two or more such lines of position can 
be used to obtain a fix - the location of the receiver. 

Various hyperbolic systems with origins in the late 1930's have 
evolved into today's alternatives. Of these, it may be said, the one 
which meets the most general operational requirements is the Loran-C 
system. The advantages of Loran-C can, in short, be attributed to 
the following. 

1. Use of a frequency band (90 to 110 kHz) which provides 
for extended range and exhibits excellent propagation 
properties, 

2. Use of a pulsed signal format, This increases average 


transmitted power and permits signal coding to facilitate 





E122 
transmitting station identification and to reduce the effects 
of interference, 

3. improved measurement precision obtained by coarse measure- 
ment of the envelope followed by fine tune measurements based 
on the signal rf 6501: 

In Loran-C, the hyperbolic navigation concept is realized in a 
somewhat complex manner. There are a number of modifications to the 
simple scheme as previously described, While the modifications make the 
system appear involved, they are seen to be well founded and straight- 
forward when systematically examined. 

A first consideration is the operation of the transmitting stations. 
These are grouped into chains whose configurations, subject to practical 
constraints, provide geometrically optimal accuracy and coverage. An 
example of a simple chain to analyze is the so-called triad which 
consists of stations denoted master, slave X and slave Y. The trans- 
missions consist of groups of pulses on a 100 kHz carrier frequency: 
the master transmits nine pulses in a group while the slaves each 
transmit eight pulses in a group. In practice, the stations do not 
transmit simultaneously. The slave stations incorporate known trans- 
mission delays relative to the master so that, within the service 
area of a triad, reception is typically as depicted in figure 1-1, 
Transmission tining is controlled so that the groups are received in the 
order and intervals shown, The transmission delays guarantee that there 
will be no overlap of the groups. In this example, the transmission 
sequence is repeated every 50,000 microseconds - the chain pulse 


group repltition interval, Different chains are distinguished by 
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Figure 1-1  Loran-C Chain Pulse Transmission Scheme 
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different group repitition intervals. Consequently, by means of a 
correlation type operation applied over a number of repitition intervals, 
& receiver can discriminate between signals transmitted from adjacent 
chains, 

The stations also employ a pulse phase code in the transmission. 
This allows enhanced discrimination between signals from competing 
chains and overcomes some of the problems introduced by skywave 
contamination. 

The term skywave refers to the multipath propagation phenomenon 
depicted in figure 1-2. The high accuracy and coverage advantages of 
Loran-C depend on the use of measurements obtained from signals 
arriving via the stable groundwave propagation mode at extended ranges. 
The reception of skywaves, whose arrival times are highly dependent 
upon unpredictable ionospheric properties, causes severe signal 


processing problems. The nature of the important properties of the 
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Figure 1-2 Loran Propagation Paths 


multipath medium, at Loran-C frequencies, is illustrated in figures 
1-3 and 1-4, It can be n that the times of arrival can be as short 
as about 35 microseconds after the arrival of the groundwave for strong 
first hop skywaves at extended ranges. Additionally, arrival times can 
be delayed as long as about 1500 microseconds for higher order hops 
under some conditions. An explanation of how these effects are 
overcome requires a further description of the transmitted signal 
format and a brief description of receiver detection methods. 

All transmitting stations attempt to generate pulses of the 
same corn which can be described as t^ a with a@ such that a peak 
occurs 65 microseconds after the beginning of the pulse. Such a pulse 
has an effective duration of 250 to 300 microseconds. An example 


which combines the two most deleterious effects of the multipath 


contamination is illustrated in figure 1-5. For clarity, the magnitudes 
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Figure 1-3 Skywave Arrival Times 








E 
NS 
2 120 Inverse 
E. Distance 
Q 
^2 100 
2 Groundwave 
S Sea Water Path 
ti e Sp : 
is St 
n 1 ' Hop SKywave 
E. 60 
un 
с 
45 Maximum Median „nd 
2 40 Noise Level == 7 M 
"d 
Helle -— A 
m 20 Minimum Median 
Er, z 
Noise level a 2 
2 “op „th Hoß \ 


Ди тл 000, 


Distance in Nautical Miles 


Figure 1-4 Skywave Magnitudes 
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Figure 1-5 Composite Received Signal 


of the separate components of the signal have been superimposed, 

The effects of the contamination can be grouped into two 
categories, The groundwaves are overlapped by skywaves associated with 
the same transmitted pulse, These skywaves correspond to low order 
hops and may be many times larger than the groundwave, The groundwaves 
are also overlapped by much delayed, high order hops associated with 
the preceeding transmitted pulse, These skywaves will be fairly 
well attenuated and smaller than the groundwaves they overlap. 

The problems of long delayed skywave interference is easily 
overcome by use of a pulse phase coding. The phase of the trans- 
mitted signal carrier is systematically varied in 180° steps from 


pulse to pulse in accordance with the scheme shown in figure 1-6, 
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Figure 1-6 Pulse Phase Code 
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As an example of the use of the code McConsideritheleffect of the 
sixth hop skywave in figure 1-5 which will cause problems by over- 
lapping the next groundwave, It can be verified from figure 1-6 that 
the code can be represented as a signal, s(t), with the property that 
s(t) and s(t-T) are orthogonal when averaged over two pulse group 
repitition intervals, where T is a multiple of 1000 microseconds, 
Since the effect of the sixth hop skywave is to create a signal 
s(t) +A s(t-T), if the received signal is correlated with s(t), this 
type of interference is removed, : 

The effects of the interference caused by low order skywaves 
are not so easily overcome, Use is made of the fact that the 
groundwave travels the most direct path and is therefore guaranteed 
to arrive before any of the skywaves associated with the same trans- 
mitted pulse, 5 ihe first portion of the groundwave is 
free of contamination of this type and can be used to make the time 
measurements, This portion, however, can not be assumed, a priori, 
to be much more than the first 30 microseconds of the pulse, 


The measurement technique proceeds in two steps, First it is 
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Figure 1-7 Envelope Deriver and RF Samplers 
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attempted to identify some point on the leading edge of the pulse. 
Since this point must be within the first 30 microseconds after the 
Defpinmane Oleic риТес ии" 15 СЈ1еах that this! s agsevere constraint 
on a pulse which has an effective bandwidth of only 10 kHz. 

A typical technique employed is to process the pulse envelope 
so as to produce a zero crossing at the 25 microsecond point. This 
is accomplished by delaying the receive em 5 microseconds and 
scaling it to produce what is represented by the dashed line in 
figure 1-7(a). This is then subtracted from the original signal - 
the solid line in figure 1-7(a) to produce the derived envelope 
of figure 1-7(b). The zero crossing is then detected and sampling 
strobes are centered at the estimated 25 microsecond point of the 
rf signal as shown in figure 1-7(c). 

With the strobes located as shown, the middle one will obtain 
a zero sample value while the other two will obtain values with a 
known ratio. If the strobes were not centered properly, different 
values would be obtained, producing an error signal which can be 
used to reposition the strobes, 

In practice, of course, the signals are noisy so that both 
Steps in the procedure must be accomplished via measurements on a 
number of pulses - as many as 64 or 128 or more - depending on the 
signal to noise ratio. Typically, the most difficult step in the 
procedure is the envelope zero crossing deriving, The zero crossing 
measurement must be accurate to within 5 microseconds or the fine 
tune measurement that follows will be made on the wrong cycle of 


the carrier. In this regard, a final characteristic of the signals 
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should be presented, 

The propagation paths involved do constitute dispersive media 
so that the signals received are not quite identical to those trans- 
mitted. The distortion in the shape of the envelope is slight - so 
much so that even though there are different propagation paths 
involved, the skywaves may be modelled as perfect echoes of the 
groundwave, The major effect of the propagation paths is the intro- 
duction of a phase shift of the pulse envelope relative to the carrier, 

In effect, the envelope and carrier propagate at different velocities, 
since the final time difference measurement is based on the carrier 
times, the carrier propagation times are calibrated for measuring 
distances, While the phase shift, for the groundwave, is typically 
not enough, by itself, to cause a cycle ambiguity in the detection, 
it does require that the envelope deriver have accuracy somewhat 
better than 5 microseconds. 

With the preceeding as background, some aspects of system per- 
es can be presented,  Inaccuracies can be introduced by any 
number of sources in the system. These can be loosely grouped into 
a few categories and statistically described, Receiver error will 
have a standard deviation of from 20 to 100 nanoseconds depending on 
the degree of sophistication, Chain synchronization includes the 
results of many effects and may be described as introducing an 
error with a standard deviation of about 50 nanoseconds, Noise 
induced errors can have a standard deviation of from about 20 nano- 
seconds for 0 db SNR to about 250 nanoseconds for -20 db SNR. The 


implications of these times can be seen from an example. If a total 
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error of about 100 nanoseconds is present, this corresponds, via a 
speed of propagation conversion,to a distance error of about 30 
meterse 
The numbers presented above for the error statistics were 
purposely given without mention of important descriptive parameters. 
For example, a figure for the error caused by atmospheric noise is 
completely meaningless without a detailed description of the 
particular receiver detection procedure - particularly the sample 
averaging or integration time involved (50 seconds in the case of the 
above figures). The geometry of the chain and the receiver position 
must also be specified: the fix obtained from lines of position inter- 
secting at right angles can be expected to have an error almost an 
order of magnitude smaller than if the lines were nearly asymptotic 
and at extreme ranges. Without going into great detail, the per- 
formance can be Simply summarized, 
1, Position errors vary from a few tens of meters to not 
more than about 200 meters within the coverage regione 
Zə The only hyperbolic navigation systems that are competi- 
tive in performance have a coverage range limit of less 
than 200 miles,  Loran-C has a nominal coverage range 
of about 1000 miles. 
3. The receiver generally contributes the smallest 
percentage of the total error. / 
In spite of the third property, improved receiver performance is 
constantly strived for since it indirectly effects many of the 


other sources of error. Chain synchronization is a single example of 
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such à source, Transmission timing is controlled by passive 
Cesium Beam frequency standards located at each station. As precise 
as these may be, there will be, over a period of time, relative 
phase drifts between the stations of a chain. These phase drifts 
are observed by monitoring stations which use the most sophisticated 
receivers and periodic corrections are ordered to maintain syn- 
chronization to within 100 nanoseconds. Clearly, if more sensitive 
receivers were available, this monitoring process could be made 
more responsive to oscillator phase discrepancies. 

Speculation on means of achieving improved processing schemes 
is tempting and many faceted. One aspect to be considered is that 
an optimum receiver need not totally restrict the measurements to be 
made in the portion of the groundwave that is perfectly free from 
skywave contamination. If, for example, the sampling strobes are 
centered at the 30 microsecond point rather than at the 25 micro- 
second point, the signal to atmospheric noise ratio would be im- 
proved. In such a case, it might result that the third strobe 
would be sampling a signal contaminated by a skywave. In many cases 
however, the skywave would be of such a low power level at the 
sampling point that the increase in total noise power - atmospheric 
noise plus skywave interference - is more than offset by the in- 
crease in Signal power. The implementation of the resulting 
adaptive time filter would require the knowledge of skywave relative 
magnitude and arrival time so that a decision could be made whether 
or not to re-position the sampling points. | 


A less optimal, but still improved, receiver would make 





— 

measurements on the skywave free portion of the signal - but at 
later times in the many cases in which skywave arrival times are 
much more than 30 microseconds after the beginning of the groundwave. 
In this case, the receiver requires only the knowledge of the 
skywave arrival time. More ambitious attempts have been made to 
implement a total power receiver. In this receiver, an attempt is made 
to digitally implement a scheme to measure skywave arrival times, 
magnitudes and phases, and then to subtract them from the received 
Signal thereby greatly lncreasing the signal to noise ratio. 

Another improvement that can be thought of would be allowed by 
a receiver capable of making precise measurements on the envelope of 
the groundwave. This would overcome the problems introduced by the 
different propagation speeds which might cause cycle ambiguity in the 
measurement, Even ifa measurement based on the envelope might 
not have quite the precision of the rf measurement and therefore 
not be of use in a navigation receiver, it vould be usefulin 
calibrating regions in vhich the envelope - carrier phase shift 
is significant, 

The investigation of an alternate method of accomplishing some 
of these goals is the subject of this report, The method - homo- 


morphic deconvolution - will now be presentec, 


1.2 Homomorphic Deconvolution 


In 1965, A.V. Oppenheim presented a theory of generalized super- 
position which could be applied to certain classes of non-linear 


e À simplified description of the theory can te given via 
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an example, Consider processing a signal consisting of two or more 
components which have been combined in some non-linear fashion. By 
means of certain elementary operations, it may be possible to transform 
the non-linear combining operation into a linear one. If the components, 
when so transformed, show a certain distinctiveness, it may be possible 
to process them by conventional linear techniques. One of the components 
can be enhanced relative to the others, or the components can be sepa- 
rated, If the inverse transformation can be realized, the recovered 
signal has been obtained by what is referred to as generalized linear 
filtering. 

The theory has been implemented with some success in what may be 
called multiplicative systems. If x, the signal to be processed, 
consists of components X4 and X, combined as in: x = əə An then the 
transformation employed might be the logarithm function (real or complex), 
What then results is, 

log x = a-log X4 t b'log xə 
so that the modified components, log X4 and log X5, have been combined 
in the ordinary linear sense. Whether or not anything has been gained, 
as well as the form of the processing to follow, depends upon certain 


properies of X4 and x_. Another class of systems to which the theory 


2 


has been applied may be called convolutional, 
If x, the signal to be processed, consists of components X4 and 


x, which are combined by a convolution, denoted x = x, x Xə, the 


2 — 
generalized linear filtering technique may be effective. The trans- 
formation might consist of, first, the evaluation of the Fourier 


Transform since it maps a convolution to à multiplication: 





7: 
Fİxİ - F İx,İ -r EA 


As before, this can now be followed by a logarithm function evaluation, 


Here, the logarithm must, in general, be complex. What results is, 


log F Ix | = log F i^] LORI EN 


Agaln, whether or not anything has been gained depends on the 


properties of x, and x,. Schafer 2 has examined a specific class 


of signals in which the X4 and x, exhibit favorable properties. His 


2 
presentation is closely followed herein. The signals may be described 


as echoed, i.e., of the form 
K 

0.1) 75 > a, s(t-T, ) 
i=0 


Such a signal consists of a basic waveform with scaled, delayed versions 
of the same waveform superimposed, This is a reasonably approximate 
description of multipath signals which can arise in many applications. 


The components X4 and x, are identified by re-writing equation (1.1) 


Z 


as a convolution of the basic waveform and an impulse train, 


x(t) = x, (+t) x x,(t) 
(1.2) x,(t) 7 
k 
x (t) E ES 24 o(t-T,) 


where ó6(*) represents the Dirac delta function 


The important aspect of the identification made in equation (1.2) 





LO TE 
is that log F EN is a slowly varying function of frequency while 
log F ES is, in many cases, a rapidly varying function. Consequently, 
tae inverse Fourier Transforms of log E E and log JEN are, in 
many cases, somewhat disjoint. A more important property is that 
- E F x is also an impulse train, Hence, by use of the 


operation ane 


| 108 F Lx] | , the convolved components are transformed 
to added components which exhibit Propert that may make separation 
possible. 

A similar operation was described by Bogert, Healy and Tukey 2 
in which the power spectrum of the log of the power spectrum of the 


convolved signals was used, The result was called the cepstrum of 


the signal, denoted x, 


nee | a Ir Bi 

Many properties of both operations are similar, In particular, 
the impulse train is mapped to another impulse train in both methods. 
A major difference is that, with the cepstrum, the loss of the phase 
information precludes any inverse transformation, Although filtering 
and recovery are not possible, some parameters of the signal are 
measurable in the cepstrum so that the technique finds many appli- 
cations, Due to the similarities involved, the transformation which 
retains the phase information has come to be known as the complex 
cepstrum while the other is sometimes referred to as the power 
cepstrum to maintain the distinction. 

On the subject of terminology; it becomes awkward to speak of 


"rapidly varying functions of frequency," as above, when the 





2 
temptation is to say "high frequency components - in the frequency 
domain." In this regard, some phrases have been developed in the 
literature of this subject and will be used to avoid confusion. In 
this thesis, the term power cepstrum will be used to refer to the 
transformation of Bogert et al. The transformation developed by 
Schafer will be referred to as the complex cepstrum, or as simply 
the cepstrum. The cepstrum will be considered to be a function of 
quefrency. Thus, if the complex logarithm of the Fourier Transform 
of a signal has a component which is a rapidly varying function of 
frequency, it has a high quefrency component. Filtering operations 
on the cepstrum will be referred to as short pass (analogous to 
low pass but in the quefrency domain), long pass (analogous to high 
pass), and comb (the usual - but in the quefrency domain) operations. 

With the advances in technology leading to high speed digital 
computation and fast A/D conversion, and, with the development of 
the Cooley-Tukey algorithm, realization of the various mathematical 
operations called for by the technique has become increasingly 
practical. Since the implementation is via digital computer, the 
theory has been developed in terms of sampled data versions of the 
convolved signals. An in-depth presentation of the theory is 
contained in Schafer * and only the relevant portions will be related 
herein, A brief overview of transform theory for sampled data signals 
is called for before the specifics are presented. 

À continuous - time signal, f(t), may be sampled every DT time 
units to produce a sequence, f(n) = f(n-DT). Such a sequence has a 


4 - Transform defined as, 





F(z) = E f(nz^ 


y =-20 
The inverse operation is accomplished via: 
- n-1 
ul S 213 9, F(z) 2 д2 


where C is some closed contour in the 


region of conVergence of F(z) 


By this operation, f(n), a function of a discrete variable is 
transformed into a function of a continuous complex variable. The 
Fourier Transform of a discrete sequence is defined as the Z-Transform 


evaluated around the unit circle in the 2 - plane, 


e со . 
F(eJ") = y fm) e” 
` y =-20 


The inverse operation is, 
1 jw jun 
fin eH Ze F(eJ") eJ" aw 
-TI 


Both the Z-Transform and the Fourier Transform vill be referred 
to in the development of the report since they allov easily accomplished 
evaluations of results, Neither, however, can be directly realized 
by a digital computer since they are functions of continuous variables. 
What is utilized is the so-called Discrete Fourier Transform (DFT) 
which is the Z-Transform evaluated at equally spaced points around 
the unit circle in the z - plane, 
N-1 
2 


ки 
F(k) = fín) e E 


n=0 
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and the associated inverse is 


N-1 jarkn 
ae = ~ E Fk) e N 


All of these transforms may suffer from shortcomings referred to 
as aliasing and leakage, Additionally, the DFT is accompanied by 
the so-called picket fence effect. These, effects and means of 
counter-acting them are developed by Ccoley et al 6 and will be 
discussed in this report only when they present problems. The DFT 
is evaluated via the Cooley-Tukey algorithm, or the Fast Fourier 
Transform (FFT). 

In the case of interest, the Z-Transform of the sampled composite 
signal is, 


X(z) =" X,(2) +X, (ә) 
The logarithm is taken to obtain, 
X(z) = log X(z) = log X, (z) + log X,(z) 
The complex cepstrum, X(n), is to be, 


x(n) = “gü köl = ES $. X(2) o 
Here the complex logarithm presents a problem since it branches 
at values of z for which ARG X(z) = d". When this happens, X(z) 
has a discontinuous imaginary part and is therefore not analytic. 
Since the cepstrum will be computed via the DFT, it is required that 
ARG X(z) be continuous for values of z along the unit circle. The 


problem is overcome via a so-called phase unwrapping technique which 
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identifies the proper branch of the arctangent function defining 
the phase, This technique is described in Appendix B. 

More specifics of the computation will be presented in subsequent 
chapters. The important aspects de the technique can be summarized: 

1. Convolutions are mapped to additions. 

2. The mapping is continuous; the inverse mapping exists 
and is continuous - hence the term homomorphic 
deconvolution, 

3. Deconvolution is accomplished via a subtraction (if 
it is known what to subtract). 

kb. The components of the original signal, in some cases, 
become somewhat distinctive in the complex cepstrum 
so that the may be identified and separated. 

The technique has found many applications since it can deconvolve 
signals when the components are unknown. It is merely necessary that the 
signal be reasonably echoed in nature with reasonably spaced echo 
arrival times. The basic waveform, either as it is or with minor 
processing, will have a low quefrency nature and decrease in 
magnitude at least as fast as ifn in the cepstrum. The impulse train 
can be made to have its louest quefrency component displaced from the 
zero quefrency line by a value corresponding to the first echo arrival 
time. In cases of large echo times, the echoes can be removed by 
a simple short pass or more involved comb filter, Conversely, the 
basic waveform can be removed and the impulse train recovered by 
use.of a long pass filter. This approach is applied in cases in 


which the multipath response, x,, is not strictly impulsive and it 
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is desired to see precisely what it is. In such a fashion, properties 
of the paths involved can be measured. 

Some applications of the technique have been to speech proces- 
sing by Schafer, Oppenheim and SCIAT s to seismology, where it is 
desired to recover and measure X5 by an. and to processing of 
measured brain waves evoked by visual stimuli by Senmoto and antics”. 

In this last application, the data is available via an electro- 
encephalogram, There is some degree of control in that the echoed 
Signal measured represents the composite response of the brain to a 
pulse-like visual stimulus. The stimulus can be re-applied, at 
fixed intervals, over a pulse repitition interval. A somewhat 
involved technique is used to average the responses obtained over 
a repitition interval so that improved signal to noise ratio is 
obtained before computing Re cepstrum. 

The similarities between this application and the Тогап-С 
case lead to speculation on signal to noise ratio improvement 
which will be discussed in Chapter 6. Considerations such as the 
effects of interference, atmospheric noise, and signal stability 
must be taken into account, Before these are considered however, 


more specifics of the problem must be developed. 
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Chapter 2 


PROBLEM FORMULATION 


In order to apply the homomorphic deconvolution technique, a 
model of the Loran-C signal as described in Chapter 1 must be 
developed. From this model the signal processing tasks can be identi- 
fied so that an appropriate solution scheme results. To this end, 
the waveform to be processed will be considered to be obtained 
from one pulse duration of the received signal, represented as 


follows: 


k - 
(2.1) R(t) öz n(t) (ox 1 s. (t) sin |, Co) -94 | 


1=0 


In equation (2.1), the zero subscripts identify parameters to 
be associated with the received groundwave, The subscripts 1 through 
k identify parameters of the first through ә” skywaves кк 
557 denotes the normalized pulse envelope of —— component 
of the composite signal. The T,'s and A, "s refer to the envelope 
arrival times and magnitudes, respectively. A represents the 
phase shift between the envelope and the carrier of the jo 
component. 

With a known carrier frequency, an equivalent method of speci- 


fying the signal is by use of the complex envelope representation 


as in equation (2.2). 


— 


k 
(2.2) f(t) = n(t) + E A, exp E İen + AS 
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Figure 2-1 Generation of Complex Envelope 


À common method of describing how the complex envelope is 
obtained is shown in figure 2-1. An analysis of the signal flow 


in the figure shows that, 


k 
f(t) = n(t) cos o A s. (t) “IL + A 
i=0 
and, 
k z 
£ (t) zu UE sd Wet + = A, s, (t) cos iT. t Al 


If the following definitions are made: 


n(t) 


n(t) sin u = n(t) cos wot 


and 


F(t) 


it follows that, 


COLE f 00) 


k 


(2.3) F(t) n(t) + da A, exp f oT, + ај) s, (t) 


li 


un 





E 
which is identical to equation (2.2). Clearly, the complex envelope 
representation for continuous signals is no more than a mathematicai 
convenience, In a sampled data form however, it can be actually 
implemented and provides a compact means of handling the quadrature 
components of the baseband signal. The realization method is as 


depicted in figure 2-2. 
Low Pass Sampler f bT) = Ер (о) 
Filter (DT sec) 
—. 


! 
| 
bol c 
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| ES 
“Sampler | er | f,(n.pr)_ fin) 
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Figure 2-2 Realization of Sampled Complex Envelope 


The component 50 is sampled and stored in the real part of 
the sample vord vhile £ (t) is sampled, complemented, and stored in 
the imaginary part. What results is a sequence of the following 


form, 


f(n) = n(n) + x A, exp 3 sen + || s. (n) 


For the remainder of the presentation in this chapter, ihe signal 
Will be considered to be noise free. The theory of noisy signals 


is developed in Chapter 5. Consequently, the signal of interest is, 
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k [ 7 
f(t) E ü A, di) J ən zh öl 5 (5) 
For notational convenience, and without loss of generality, the 
magnitudes will be normalized with respect to Лу and the phase 


notation will be simplified: 


A, 
ДИ = = ? ј = О, 1. 5. 
О 
2, = wor. 4+ o; 
B. = A; e1; 
di l 
Hence, 
k R 
(2.4) f(t} = = 3. s.(t) 


In Chapter 4, distorted signals will be considered. Until then, 
it will be assumed that s, (+) = s(t-T,). In this case, equation 


(2.4) becomes 


k 
f(t) = 2 B s(t-T.) 
i-0 y 


This is the form of signals for which the homomorphic deconvolu- 
tion technique has been developed. Additionally, it is possible to 
identify the tasks to be performed at this point. As presented in 
Chapter 1, it is desirable to measure Bos the groundwave parameter 


from which time measurements are made, and T., to measure the difference 


o? 


in propagation speeds of the envelope and carrier. Additionally, 








—” 
measurements of Tn and B4 would be helpful. Finally, measurements 


Df all f, 's, cand İB 


i NE would be desirable since a total 


power receiver might then be feasible. 

With the problem thus defined, solution techniques can be 
considered. To begin, consider the following to be a reasonable 
approximation to the received Loran-C pulse envelope: 


s(t) - 5. u_,(t) ; co Gay x 10° or 


where u_, (+) is the Heaviside unit step function 


Initially, make the following assumptions: 


J. = 0 

T. = n, DT 

m an integer 

DT = the sampling time 


The sampled waveform in this case is: 
k k 
on) ee x B. s(n-n, ) = s(n) @ 2 В, ö(n-n,) 
i-0 i-0 


where (9 denotes a discrete convolution and ö(n-n, ) is a unit sample: 


i n=n 
6(n-n, ) = 
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Therefore, 
f(n) = s(n) ® p(n) 
where 
k 
eG) = “əə B, ö (n-n, ) 


In the case of k = 1, i.e., in the case of a signal consisting of 


the groundwave and one skywave, 


p(n) = Sh)+ B, 6 (n-n, ) 


163 [B4 < 1, the complex cepstrum of p(n) is (see Appendix A for the 


details of the computation): 


m 


oo B 
$() - x (077 —i £(n-nmn,) 
m=1 2 
whereas the cepstrum of s(n) is, 
= + (-1 mH Am 
2.5) $(n) - -x6(n) + = 5 e 6 (n-m) 
=1 


Since the complex cepstrum operation has the property of mapping 
convolutions in the time domain to additions in the quefrency domain, 
the composite cepstrum consists of S(n) + D(n). In this case, a 
simple estimator can be made to Scan the cepstrum sequence and 
determine n4 from the location of the components of the impulse train. 
À simple comb filter could then be employed to remove all traces of 
the echo from the cepstrum - only slightly effecting the component 
due to s(n). Before filtering p(n), an estimate of EM and 2, 
could be made from the components of the impulse train. In such a 
manner, estimates of T,, |B, | cud A are obtained and all skywave 


contamination is removed, 
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Vhen T, is not an integer number of sampling times greater than 


T., P(n) no longer has such a simple form. The components of the 
sin” 
x 





impulse train become terms which spread throughout the cepstrum 


sequence aS in equation (2.6) (See Appendix A). 


00 ph a ruri +A y 
A m+1 1 ləri 
(2.6) p(n) = e (-1) m min-m n, +A; j 


an integer 


( 


where ny 


and uğu DT-(n, +A,) 
- 0.5 <A, = О. 

In this case the filtering becomes less straightforward. A simple 
comb filter can not be effective. In the distortion free, noiseless 
case however, the cepstrum consists only of S(n), which is known, 
and p(n), which is known to within three parameters: Bi, nq, and A,, 
so that an estimator can be developed to identify p(n). When this 
identification is made, the deconvolution can be accomplished via 
a’ subtraction. The form of the estimator is developed in Chapter 3 
and the effects of signal distortion and noise on its performance 
are presented in Chapters 4 and 5. 

The problem is further complicated when To is not an integer 
multiple of the sampling time. In this case, the s(n) Of equation 


(2.5) has another term added to it, (Appendix A) 


S(n) + CES — ; nf£O 
^ 


(2.7) s'(n) = 


s(n) ; n20 
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The presence of this new term makes the estimation procedure described 
above suspect. Now, there are two unknown components in the 
computed cepstrum. Since the forms of the components are known 
however, a modified estimator sa still be effective, 

A summary of the problem as formulated thus far can be made 
via figure 2-3. The received signal is processed to produce 
quadrature baseband components, These are sampled and digitally 
combined to produce the sampled complex envelope. The complex 
cepstrum is computed, From the cepstrum, signal parameters are estimated. 
One of the uses of the estimated parameters is shown: the groundwave 
phase estimate provides feedback to the local oscillator to permit 
phase tracking. An important aspect of this particular use is that 
it provides a method of comparing the performance of the homomorphic 


filter with other receiver techniques, 
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Chapter 2 


DECONVOLUTION OF NOISELESS SIGNALS 


As developed in Chapter 2, it is desired to estimate waveform 
parameters from the cepstrum of the received signal so that the 
separation, or deconvolution, can be achieved. In the distortion 
free, noiseless case, the technique is as presented in this Chapter. 

It is presumed that, with one skywave, the received signal will 


have a cepstrum of the following form: 


f(n) = $'(n) + (n) 


7.1) 


s(n) + @(n) + 5(n) 
where, from equations (2.5), (2.6), and (2.7), 


00 m+i 
(a) 8(n) = -a 6(n) + > Bt C1) T АК 
n= 


O 


З 
| 
о 


7.2) (b) gün) = ayid onun 


sin nİn-n(n, “4, ) 
1i |n-T na tA 


00 B 
6757” 





In the estimation procedure it is desired to use the knovledge of 
the form of the above sequences in processing the received signal. 
Consequently, it will be necessary in the development to refer to the 


above sequences and corresponding components of the received signal 
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simultaneously. For clarity of presentation therefore, the notations 
s(n), g(n), p(n), and f(n) will be used to refer to the predicted, or 
ideal, sequences as defined above whereas the received signal will 


be described as, 


r(n) 


İl 


s.(n) 6 e(n) @ p (n) 
(3.3) 


a) - $.) + En) + $,(n) 


The subscript a signifies actual components of the signal as 
received, In the distortion free case, the only difference between 
the components of (n) in equation (3.2) and f(n) in equation (3.3) 
is due to computation error. 

A third set of sequences will also be referred to. Estimates of 
signal parameters Ao, Ai, nə, and By will be made from the received 
Signal. From these, estimated sequences will be generated, The 
estimated sequences will have the form of those in equation (3.2), 
but with parameters as estimated from the sequences: of equation (3.3). 


Thus, the estimated sequence of £(n) is: 


oJ 
^ A 70 
A A n 
g(n) = g(n) A =A = 
О О 0 ; п = О 
where Ao is the estimate of A, obtained from rM 


dei Sequence of Estimation 


The received signal consists of the terms 8. (о), & (n), and p, (n) 
Which are combined to produce the cepstrum sequence shown in figure 


3-1. 








Figure 3-1 Cepstrum of Simulated Received Signal 
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It can be seen that there is insufficient seperation between the 

sequences to allow for precise estimation of the various parameters 

directly. Consequently, the estimation is accomplished in a number 

of steps; systematically estimating the parameters one at a time. 


The procedure used herein is to make the estimates in the following 


order. 
IS 
1. S(n) - an estimate of the received signal basic 
waveform, 
zə (n) - an estimate of the small linear phase term 


not removed before computing the cepstrum 


p. 
3. Pin) - E 
(a) n, - a coarse measurement of the 


skywave arrival time 
(b) 0m - a precise measurement of the 
skywave arrival time 


pn 
(c) B, - skywave magnitude and phase 


The procedure begins by utilizing the fact that a good a priori 
estimate of S (n) is available, to wit: S(n) of equation (3.2). 
Therefore the first estimate is Ze) = §(n). This estimated 
sequence can be subtracted from the received cepstrum to generate a 


new sequence denoted x(n), 
^ ^ A ^ ^ : A 
AO AS s(n) = ео A) € (n) 


where 
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In this chapter İt can be assumed that this error sequence, nm 
is negliglbly small, The procedure as thus far discussed was accomplished 
on à simulated received signal and the resulting X(n) is shown in 


figure 3-2. 





Figure 3-2 Simulated X (n) Sequence 


It can be seen that there is improvement over the situation shown in 
figure 3-1, but that the components of 8, (n) and p, (n) still interfere 
so that the desired estimates are still not readily obtained, More 


precisely, it can be seen that for quefrency values greater than 





terms of p, (n) interfere with one another 


(n,+A,), the various n a 


to some extent, For positive quefrencies less than (n,+A,), there is 





ES 


SIN x 





interference between & (n) and the first term of p, (n). 

What makes the estimation possible is the fact that for small 
negative quefrencies, p, (n) is reduced enough so that it only slightly 
interferes with & (n). The next ep in the procedure therefore 


will be to estimate A) from the negative quefrency portion of the 


cepstrum. 
5722 Estimation of 8, (n) 
Since &(n) = —. ES 150 О, 


it follows that 


Ao di 8, (-1) 
дә 
J = E &, 72) 
I 
50 . 2 (-3) 


иә 
where Ay Will be the estimate of Ao: As described above, it is true 


that, 


$(1) « & C1) + $ (-1) + Ê(-1)= £(41) 


r(-2) 


g (-2) + $ (-2) + €(-2) = &(-2) 


since the B, (n) sequence is small for this range of quefrencies. The 


estimate of Ay can therefore be obtained fron, 
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1 + -5 Ea m 


^) 
With A, so estimated, the sequence &(n) can be generated and 


subtracted from X(n) to produce the new sequence y(n): 


Hn) = &n) - An) 


II 


where 


(n) 


It 


Sie Be o 


n 


The procedure was applied to the X(n) in figure 3-2 to produce 


the y(n) of figure 3-3. 
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Figure 3-3 Simulated y(n) Sequence 





ə 
At this point in the procedure, since Ty = DT (ng*À9) , the 


estimate of To can be made: 


/ə AZ zə 
To = DT (n¿+4,) 


where ffi. is deduced in the phase unwrapping computation as described 
in Appendix C, 
The reason for the particular estimator used in equation (3.4) 


can be seen from figure 3-4, 
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p(n) 
envelope 
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Figure 3-4 Comparison of &(n) and p(n) 
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Both sequences & (n) and p, (n) decrease essentially as fast as 1/n 
for negative values of n. Near the zero quefrency value, the & (n) 
term clearly dominates the composite signal. As n becomes more 
negative however, the distinction between the two components is not 
as great, Consequently, only a few samples should be used to make 
the estimate and even among these few, the samples obtained from the 
more negative quefrencies should be less heavily weighted, This is 
precisely what is accomplished in equation (3.4). 

The next step in the procedure is to estimate the parameters 
of p(n), i.e,, A,, ny» and B4» from the sequence Y(n). This is 


accomplished in three steps, beginning with the estimation of Nye 


3,3 Estimation of р, (п) 


The values of y(n) are examined for quefrencies corresponding 
to times betveen about 30 microsecönds and 60 microseconds, The 
largest value of ¥(n) in this range of quefrencies will be assumed to 


sin = term in 





correspond to a sample in tne main lobe of the first 
D (n). Consequently, the value of n at which this maximum occurs 


Will be the estimate of n The situation would typically be as 


Jt 
depicted in figure 3-5, Since the values of Y(n) will, in general, 
be complex, the plot in the figure corresponds to either the real 


or immaginary part of ¥(n). 
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Figure 3-5 Expected Form of pa(n) 


Once n, is determined, B4 and A, can be estimated as follows. 
Consider the values of y(n) in the vicinity of n= n, to be due 


öniy to tne first 





= X function in the expansion of $(n). As an 


example of the validity of this assumption, consider the values of 





the first and second terms in the sum of equation (3.2b) for n = ny: 
ş y 2 9 a 
A E". Sin "ny - (ny +14 )) i B, sin ија, 2(n, +A, ) 
T e 9 9 
: 2 
p B, sin TA, | B, sin T (n, - 2h, ) F 
A, 2T - n,-2A, 


The ratio of the magnitude of the second term to the magnitude of the 


first term will be, 





5 





A: DM |sin 214, | Ay sin A, 
2 [ny 72A. | es ¡sin TAy | = n, - am, | Pa 2 sin TA, 





Exponentlal weighting , as described in Appendix A, ensures that [B is 


less than 1. Additionally, the term 





sin 2TA, 
2 sin TA, 
is less than 1 so that the ratio is less than 
|Ai 
5.5 1 
—Дһ UN E ; LOrE 1 
Ы - 2A | c E HA. 1 


since АЈ ci 


Hence, for reasonably large values of ny» the approximation 


is a good one and allows the following: 


да — B, sin ulu -2- (ns +A, )| R B4 sin T(2+A,) 
A Tin, -2- ny +A) n(2+A,) 
— B, Sin mn, -1 - (ny +h, )! By sin T(1+A, ) 
iita Tin, - 1- ny +54) T(1+A4) 
bo An.) B, sin İn, - (n+) B, sin 7A, 
| 1 "in, n4*À4)] 
“ Los B4 Sin Tn, + 1 - (n,+A,)! a B4 Sin T(1+A, ) 
1 min, + 1 - ni*À "LEA, ) 
" = B, sin "mn, +2- (n, +A, )) - B, sin 1 (24A, ) 
5. " uz +2- (ny +04) m(2+A,) 


The solution can proceed in two steps: estimate A, then Bi. To 





ups 


accomplish this, define the following variables: 


d, = y(n,-2) - A4 
¥(n,) 2 tA, 
y(n,-1) > A, 
— E 
2 $(n,) i +A, 
(3.7) A 
y(n, +4) A4 
d = 
3 (n4) 1 = бә 
y(n, +1) = A, 


It can be seen that four estimates of A, can be obtained from 


the d, "s, The estimates obtained from d, and d, however, are made from 


2 2 
values of y(n) which are closer to, or in, the main lobe of the first 


sin x 


= term, They are therefore less sensitive to the effects of 


sin x 





interference from other terms and the residuals of the 8. (n) 
estimate than the estimates obtained from ds and de It is, however, 
desirable to base the estimate on as many measurements as possible 
so that a compromise would be to make a final estimate based on a 
weighted sum of the four individual estimates, The procedure is: 
1. Compute the numbers da, d, ds and dy as in equation 
(3.7) 


2. Make four individual estimates as follows: 


A e diu 
11 5 
— _ - d, 
Aio = 1+d 
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3 

Ay, 2577: 
3 

- d 

m a a 
4 


where ди denotes the estimate of A, based on the 
measurements which generated dy" 
3. Weight these estimates by the corresponding magnitudes 
of the measurements generating them and perform an 


appropriate scaling to arrive at the final estimate: 


4 n? 
a ду | Àir 


ba 4 
“m 
EN. 
When this estimate of A, is obtained, equation (3.6) can be 
re-applied to estimate Bye To do this, make the following individual 


estimatese 


8, = TÜ 521) (1) 
sin T(1 + A,) 
/? 
C S E $ (n, ) 
13 sin m, 21 
AL 
"(1 - A) N 


Ы3 2 
I 
o 
po 
3 
EI 
— 
t 
— 
< 
— 
I 
P 
+ 
— 
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14 





These are combined to produce the final estimate: 


5 . 

uu SR B 
N M a, +3 ) | ik 
By = 7s - 

= [$0 3-11 

Kel 


With the estimates of n¿» Ay» and Bi, the estimated sequence 
5) can be generated. This a) can be subtracted from the 
received signal cepstrum, r(n), to produce 260 If the inverse 
DFT of É(n) is computed, what would result would be an estimate of 
the complex log of the groundwave, From the zero frequency point 
in this sequence can be measured 2, - the phase of the groundwave, 


If there are more than one skywaves present in the received 


sin x 





waveform, 9 (n) will have additional terms centered at 


quefrencies of (n,+A,), 2(n,+h,), e (n,+A,) 55— and 

at various combinations of (n, +A, ) and (n,+A,) as described in 
Appendix A. Once the terms due soley to the first skywave are 
removed however, the lowest quefrency Sinx term that remains is 
due only to the second skywave so that the above estimation procedure 


can be re-applied to estimate nə Az and B This can be done as 


2 
often as necessary. The entire estimator can be described as in 


Youre 3-6, 
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3.4 Performance of Estimators 


Before the results of applications of this technique are pre- 
sented, a reason for the ad hoc nature of the estimators can be 
given. Due to the sequential nature of the estimation scheme and 
the manner in which the received signal is mapped to the cepstrum, 
the estimators are called upon to act in the presence of somewhat 
Structured interference. As will be discussed, the interference 
Structure can be anticipated so that the estimators can be designed 
Eomprovide Some degree of protection agalnsıt it, To show this, it 
should not be necessary to work with the complicated estimators as 
developed. Instead, a similar, but much simpler, estimation 
Situation can be considered. This will allow a clearer presentation 
of the important property of the estimators. 


As a rough model of the sequence to which the A, and B, estimators 


1 1 


are applied, consider the sequence z(n) shown in figure 3-7, 





Figure 3-7 Simplified Signal to be Estimated 





7 
The sequence is knovn to be of the forn shovn, but the value of a 
is unknovn and to be estimated, The important manner in vhich 
z(n) is similar to D(n) is that the signs of the z(n) sequence 
values have a particular code to them: the values of z(n) are positive 
for odd negative values of n and even positive values of n. This 
code is similar to that encountered in p(n) sancı En D(n), there 
are two Samples in the main lobe of the first Ba term, causing 
a Similar sign reversal. It is desired to consider the effects of 
the presence of structured interference on the estimation of a. 
To do this, consider first, the estimator of equation (3.8) which 


is of the form used for A, and Bye In effect, the estimators are 


matched to, or instep with, the code of the z(n) sign reversals. 


a, = -2 2,(-2) 

2, = 2.5 
(3.8) 

a = = z (1) 

ay = 2 z (2) 


Now suppose that the actual sequence, z (n), consists of 
the anticipated z(n) as in figure 3-7, plus a constant interference 
sequence, w(n) = b. The sign code of wln), for n = -2, -1, 1, 2, 
iS + + + +. This can be described as being orthogonal to the sign 


code of the estimators, which is + - + =. 
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N 


MAS caen z (n) z(n) + w(n) = a(n) + b, so that, 


a, = -2(-— +b) = a-2b 
2, = (a+b) = a+b 
em T Zaf+ib) = a-b 
3 

PE 40) = a+2b 


Hence the individual estimates have complementary error. It should 
be assumed that [a | =2|Ы and a> O in making the following 


argument about the combined estimate, 





Since a »|2tj it follows that, 





Eco]? [o +] - (-» 
ə Gİ - |(a + »)| se E) 
2,0) = |e- Ы) = (а – 5) 


* 


|а, (2) x (+9 = = : 
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From this, it also follows that, 


(5 - b) (a - 2b) 4 (a s b)” 4 (a - b)” 4 (ğ 4. b) (a + 2b) 


Vi 
aat 
> – Ы) +'(a + b) + (a - b) + G + b) 
2 2 2 
(3.9) a - ри =a + 2—— 


As an example, suppose that b = a/5: 


= a +t —— = 1,08a 


The estimate for this technique is 0.08a whereas if only z (1), 
for example,were used, the error would have been b - 0.23. If b - „Ía, 
the errors are 0,02a and O,1a respectively. If it happens that the 
interference has the same code as the estimator, it can be shown 
изи тһе error is 0.36a for Ы = ,2а and 0.17a.for b z Ja 


If w(n) 2 (-1)! v», the individual estimates will be, 


ə pa = = 

ај = -2(- > dE b) a - 2b 

a = (a - b) = a-~ b 
2 

2, = S (-a - b) = atb 

N a 

aj 2 ( “nal b) = a + 2b 


Once again the errors are complementary and, it can be shown, the 
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total estimate will be as in equation (3.9). Figure 3-8 shows the 
variation of the estimate error as a function of b for an interference 
sequence of constant magnitude. The three plots are of: estimate 

from z (1) alone, four step estimate - interference in step with 


estimators, four step estimate - interference orthogonal to 


estimatore 
ə 
651, / 
94; 
va 
- 
7 
Interference P 
in step with 
4 pt estimato A 
H 
O A 
E tept A 
0 25a estimator interference 
"n 7 orthogonal to 
2 p 4 pt estimator 
(9 z 
"a 
- 
a 
PA 
24 
2, 
b 
2a > a 


Interference Level 


Figure 3-8 Estimator Performance Curve 
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In the actual case, the estimators are a bit more complicated 
than those of equation (3.7) and the symmetry is not quite complementary 
as described. The example does however, illustrate the improvement | 
obtained by centering the measurements around the main lobe of the 
sinx functions, Two cases wherein this property is helpful can 
be presented. 


The use of the estimated sequence g(n) will generate some 


Small residual: 


A Al 
en) = & (n) - gta) 


& 


En 
- LL (ap - (a) Y 


This is one source of interference when the A, and B4 estimates are made. 
Although the sequence is not of constant magnitude, the change is 

only fron (€a,)/6 to (En 7/10 in the quefrency range of interest 

for n, = 8. This is à reasonable enough approximation to a constant 
magnitude to produce some improvement. A more important case, since 

the magnitudes involved will typically be greater, is the interference 
from the second sinx term, As developed in equation (3.5), the 


magnitude of this term will be fairly well reduced for quefrencies 


near Nye The situation can be examined a little more closely to 


sin x 





see that, if the largest value of the first term is denoted a, 


where B, sin T. 


125 — 


then the interference sequence, denoted w(n), due to. the second 





— 

















ean X function has values, 
2 
B, sin vA A 
jw (n,+2)| - |2 + A : 2297” EN 2(n E: E 2) 
B cl 151 
дк 
B sin nA, A 
w(n,+1)| = l 18 1 
1 2n(n, +A, al) il 2(n, +A, 21) 
: Big sin TA, = |» | A, 
= = p 
may 2) S YA, 32) бу Да) 
If B, = 0,8, ny = 8, and A, = 4, these values range from about 


0.03la to about 0,019a. Once again, the symmetry is not perfect, 


but it is close enough so that some improvement can be gained, 


3.5 _ Application to Simulated Signals 


In the case of the generated received signal used for the 
example depicted in figures 3-1 through 3-3, the following parameters 


were used: 


DT = 5 (microseconds) 

Tə 11.5 

Т, = TQ 23995: my = 8 A, = -01 
P4 ES 0.8 

fp = 0 


When the technique developed herein was applied to this signal, the 


resulting eStimates were: 





E 


re, 
2 = 
n4 8 
A, = -0,096033 
" 
DEZ 0 O5 108 
B, = (0.803417, -0.000881) = 0.803417 - j0.000881 
Ø, = 0.000686} radians 
ni na 
The phase error in this case, [A - dol m Po: 1520.039 in 


degrees, At a carrier frequency of 100 kHz, this corresponds to a 
time error of about 1.1 nanoseconds. The procedure was repeated 


for a variety of cases and the results are presented in Table 3-1. 










, Al Ad 
Pol To | T} |B 24 error 
nanoseconds 


= 


1 
, 000004 
=. 0001 
-.000 6 
— 001 
11,5138.0 rape 37.9698 |(.803411, Eo 
.00019 
8.5 140.0 gm (2795314, 380 
00024 


HARCA 57-27: (4799695, 
,000 

8.,5[38.0 |.8] 8.3650137.93861(,808735, | -,001570 22.50 
00013) 


Table 3-1 Results of Simulation: No Distortion, One Echo, Zero Phase 
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ә 
— 
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(Note: For convenience, the T, column in the tables actually represents 


skywave arrival time relative to the groundwave; i.e., (Table T.) = Tu 





The results shown in the table reflect the performance of the 
estimators on signals which were purely real to begin with. The 
phase errors would therefore define the precision available in the 
loop shown in figure 2-3 in which the homomorphic filter-estimator 
serves as the phase sensitive detector. 

The results of the simulation for non-zero phase cases are 
shown in Table 3-2. What can be seen in these cases is that the 
B. "s are essentially -.19723 radians out of phase with the true values, 
This can be explained by examining the Fourier Transform of p(n): 

: ; B à 
le) = m «2, e 74 - m (14 5 e791) 


In other words, all magnitudes and phases in the complex cepstrum 


are measured relative to the groundwave magnitude and phase. 


8 : 


"19722 mo Е4 77: 25 1 ии m. (. 783872, » 197302 5112 


a 7 5 ES 0 bö 25 — 9292 "792381, «202399 


Table 3-2 Results of Simulation: No Distortion, One Echo, Non-zero Phase 





The cases in which the skywaves are larger than the groundwave 
are considered next. As discussed in Appendix B, exponential weighting 
must be used in these cases, The only thing that changes in the 
estimation procedure is S ihe a priori estimate, must also 


show the effects of the exponential weighting. This is easily 
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accomplished as shown in Appendix A, The results are presented in 
Table 3-3, The extra column, Bi» gives the effective skywave coefficient 


relative to the groundwave after exponential weighting. 
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2,0 
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78540 110.0 20. 775 O53317, | 7284751 C0 
eR ,000722 





Table 3-3 Results of Simulation: No Distortion, One Echo, Skywave 
larger than Groundwave 


In the cases presented thus far, either the groundwave and the 
skywave were in phase, or one of the two was at zero phase, For 
completeness, Table 3- presents the results for different, non-zero 


phases o 


.19723110.0120.0 |(2, 3) | (.39233, |10.0002 |39.9993| (.392375, „064 
“e 8828 


TT 58821 
(19723 110.0139,5 (4, -3) | (639233, | 9.9412 /39.4472 ¿201864 17.36 


| m e 8828 1 
419723 139233 9.8134 137.99941(.38774, 1.199813 14.10 
- ,58828 - , 59804 


Table 3-4 Results of Simulation: No Distortion, One Echo, Skywave 
| and Groundwave at Arbitrary Phases. 
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İn all of the above cases, the final estimate of the groundwave 
cepstrum sequence, А (о), was short pass filtered after 150 micro- 
Seconds before taking the inverse DFT to measure Bo As discussed 
in Appendix A, this is done to facilitate filtering in the multiple 
Skywave case, Table 3-5 shows the results of the estimation scheme 
when applied to a signal with two skywaves. The results are not 
quite as good as in the single skywave case. Nevertheless, it is 
felt that, if warranted, more elegant estimation algorithms can be 
developed. The point of the presentation thus far is to demonstrate 
that, in the distortion free, noiseless case, measurement precision 


on the order of a few nanoseconds is possible, The effects of 


Signal distortion can now be examined, 
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Chapter 4 


DISTORTED SIGNAIS 


The results of the preceeding chapter depend heavily on the 
assumption that the received signal is identical in form to the 


A 
sequence f(n) defined in equation (3.1): 
2 A A ^ 
(ne sm Hen) Et (n) 


In practice, it can be expected that & (n) = e(n) since this component 
is merely due to the small portion of the linear phase not removed 
before computing the cepstrum, There remain however, two sources 

of error since S (n) and p, (n) may not be of the expected form. The 
sensitivity of the estimators to these errors will be examined in 


this chapter. 
4.1 Distorted Basic Waveform 


In the estimation scheme of Chapter 3 it was assumed that a good 
a priori estimate of S (n) was available. What is meant by good might 


be that some criterion such as equation (4.1) is met. 


AJ 
S (n) - S(n) 
EEXONE E C för ail n such 


that S (n) = О 


ê (n) 


(A, 1) $ (n) 











where € is some appropriate bound. Upon close examination however, it 
can be seen that what is of importance is the effect of the residual 


^ 
Sequence, € (n), on the sequence X(n). 
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x(n) 


6(a) - $a) - $a) s £.(a) 4 Ф. (о) = 20) 


= £ (n) +9 (n) +E (n) 


Consequently, what is of concern is the ratio, 


€ (n) CN) 
(4,2) MES Ta - ken ; föorgall nol 


interest 


As an example of the distinction, the ratio of equation (4.1) 
might have a value of about 0.01 for some value of n. As this stands 
by itself, (n) might therefore be considered a reasonably good 
estimate. However, if S (n) is 100 times greater than p, (n) * & (n) 
at this value of n, the ratio of equation (4.2) would equal 1. The 
estimate might then not be considered so good - particularly if the 
value of n in question happened to be near nje | 

Conversely, if the ratio of.equation (4.1) is quite large for 
some quefrency value not considered in the estimation procedure, or, 
for some quefrency at which p, (n) + e, (n) is much larger than $ (n), 
there may be no problem, The criteria for judging the appropriateness 


of the a priori estimate would therefore be as in equation (4.3). 


A 


(4.3) C (n)«c& (n 3  n=-1,-2 -3, -h 


£ (a) << 9 (n) : nen 


22 


1 


Distortion of S (n) Will be introduced by the propagation path 


and by the receiver processing. It may be assumed that the propagation 
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distortion of the groundvave is negligible, Not so negligible 
however, would be the effects of receiver filters. In addition to 
performing the low pass filtering in the demodulation scheme shown 
in figure 2-2, the filters are called upon to suppress atmospheric 
nolse and interference from other navigation or communication systems 
operating in the low frequency portion of the spectrum. Since no 
general linear time invariant filter can be used to model the effects 
of these filters under all conditions of operation, what will be 
presented will be a method of analyzing the technique which could 
then be applied to any specific implementation. 

suppose the received signal is processed by a Single pole low 
pass filter (similar arguments can be applied for double or triple 
pole filters). A model of the basic waveform, i.e., Ei 


without echoes, to be sampled would then be, 


s (t) ECC) x h(t) 
where noo .. uq (t) 
so that S. (n) =-S(a) + hin) 


and € (n) = $ (n) - 2 = S (n) - S (n) = 34 


In this case, the cepstrum of the filter response is easy to 


compute: 


H(e") = x o... h - UG) 2 
n=0 





17 








e co e Pn 4 
log Hc 15) =i F = ey Гог В = О 
n= 
and 0 : otherwise 

^ ^ 

Ê (n) = h(n) = 

= 
oo eae 
n 
“ 4 


The fact that the cepstrum is zero for negative quefrencies will 
be true for any minimum phase filter - any linear phase term being 
removed in the phase unwrapping procedure. Since the A, estimate is 
made from the negative quefrency values of the received signal, it 
can be seen that the accuracy of the estimated sequence g_(n) will 
not be affected by the presence of the h(n) term. The values of 
€ (a) for quefrencies near n, however, will cause problems. The 
values of € (n) and B(n) in this region are given in Table 4-1 for 


comparison. The parameters of the signals in this case are: 


BOX =E 2/65 
Haa 8 

1 = -0,1 
B4 = 0,8 
DT = 5 
















€ (n) P(n) 
13657 
n,-1 .11518 ,0792 
08523 | -.0791 


Table 4-1 Comparison of A Priori 
Estimate Error and $(n). 
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ə 
€. (n) 
Ercarbs "$n is too large in this case to permit estimation 
of A, and Bye The problem can be overcome if the response of the 


filter can be presumed to be known reasonably well. The a priori 


estimate of S (n) can then be modified: 


S (n) = Ә(и) + dn) 
so that € (n) - s (n) - &(n) 


= §(n) + h(n) 2 sn) ee h 


(n) 


A X 
= h(n) - h 


(n) 


In this case, an estimate of the form 


О : otherwise 





; n»0 


Where y x p might be used, The error therefore would be, 


€ (n) — i-|. ~ zn | "B I0 


For Y = 4, with the other parameters as in the preceeding example, 


the error values would then be as shown in Table 4-2, 


¢ 





| ' 
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P(n) 


01176 


Table 4-2 Comparison of Revised 
КиһетОта Estimate 
ErrorBand pn) 









Within this range of quefrencies, € _(n) < 9 (n). Further 
notice should be taken of the fact that the € (n) terms are of 
relatively constant value and of constant sign. As developed in 
Chapter 3, the individual estimates of A, and B, will therefore 
have nearly complementary errors, Consequently, the total estimates 
of A, and B, are relatively unaffected by the incorrect a priori 
estimate of P. 

if he correct value2of Y is used, 1.8., Y. By but an incorrect 


filter form is assumed, the results are quite different, Suppose 
O 


It can be shom that the resulting cepstrun is, 


0 s. n<0 
^ 
2525 UD ; n = (О 
see 
2 п = 0 





n 





E 


Therefore, 


o 
5 
A 
o 





€. (a) ә -p | g. ПД 0 
-Bn 
x :5n- O 
mr J 


This will produce the same situation as that depicted in Table 4-1, 

In summary, it can be said that the a priori estimate of $ (n) 
must be changed to reflect the effects of the responses of any filters 
used. The estimation procedure of Chapter 3 will be relatively 
insensitive to errors in filter parameters, €.8., pole locations), 


but fairly sensitive to errors in estimated filter form, 


4.2 Distorted Echoes 


If the received signal is processed by a linear time invariant 
filter, the effects will be as described above. Such a filter will 
not change the echoed property of the signal, i.e., B (n) will not 
be changed. Distortion of p. (n) can only be introduced by Air acrenc cs 
in the propagation paths, The distortion can be modelled as follows 


in a single skywave case, 
r(t) = s (t) * B, s (t-T,) 255) 


к(е9Ҹ) - 5 (е9").с (еУҸ) mi +3, ¿wn LU, H(e JW) 
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log Hs ) log SM + log с̧ (е7") 
a .. yA E Zr 
ə m) "= “= (et) 
kei 
for B, H(e") | < 1 
A ^ ^ ^ 
r(n) - s,(n) * &(n) + p(n) 


3 k әтәк Са +A,) 
(n) = Fİ nü әт. 


A 
where p k 
k=1 


a 


The resulting 5, (n) can be anticipated by recognizing that its 


components are inverse Fourier Transforms of products. Hence, what 


sinə 





was the first term in the expansion is now convolved with 


Sime xX 





h(n) = h(t) . The second 
t-n:DT 
twice and the third term is convolved with h(n) three times, etc. 


term is convolved with h(n) 


The results of these convolutions are combined to form p, (n). In 


the cases considered thus far, h(t) has been an impulse function 





25 the sin x 


terms appeared intact. If h(t) is only approxi- 


mately impulsive, the first mcum term will be only slightly 





sin x ¿ : 
x terms however, will be convolved with 


distorted. Higher order 


h(n) à number of times and will tend to spread out - only loosely 





functions, 


| sin x 
resembling x 


In the Loran-C case, h(t) is nearly impulsive. Additionally, 
Since the cepstrum is short pass filtered, the effects of the many- 


5. terms vill 





fold convolutions vhich distort the high order 





"m 

only slightly be felt. The selection of an appropriate model for 
h(t) is confounded by the fact that it is so nearly impulsive that 
it has not been possible to measure its form. Perhaps the only 
correct model would be that the skywave and groundwave are identical 
to within 1% along the leading edges of the pulses. 

Mitholehtan appropriate form of hit) is not available) it is 
desired to test the sensitivity of the estimators to a possible 
1% distortion. What will be used in this case is a groundwave of 
the form t m and a skyvave of the form t* a Ifa = 2/65 
and B = 2470; and if the signals are scaled to achieve a value of 
1 at their respective peaks ( 2/0 and 2/8), the skywave will have a 
magnitude of about 765 of the magnitude of the groundwave at t = 65, 
(with both signals beginning in coincidence). While the h(t) that 
would produce such distortion, ices, rt (fe), ïs not 
claimed to be a particularly realistic model of what actually occurs, 
it does serve as a clearly large overbound on the amount of distortion 
to be expected. 

The case described above was simulated and the rls are 
presented in Table 4-3, While the errors are, understandably, 
higher than in the cases presented in Chapter 3, they are the results 
of distortion which is much more than can be expected to be introduced 
by the propagation path. As long as the propagation path filter is 
nearly impulsive, the effect on the cepstrum estimators is minimal. 

The results presented in this chapter will be discussed further 
in Chapter 6. It can be said however, that, with the assumptions made 


in the presentation, signal distortion only slightly affects the 





Eo 


performance of the homomorphic filter as used in this application, 


| 
iO, 


(.32977, 5—— (.31743, aa 
s - “4948 - 44-7539 
© 33627, i in .34000, 1.18663|16.9 
~, 50421 - ДӘ 
»19723110.0/38.01(2., |(.35652, | 9.8317/38.0789/(.35804, | .18418/20.8 
-2.)| -.53459) -.55027) 
Table 4-3 Results of Simulation: Distorted Echo 
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Chapter 5 


NOISY SIGNALS 


In this chapter the effects of noise on the cepstrum estimators 
will be examined, Before quantitatively presenting the effects of 
the noise, an attempt should be made to develop a theory which 
outlines the results to be anticipated. The development begins 
with an explanation of a matrix notation which seems best suited 
to a clear presentation of the elementary aspects of the problem. 

If a continuous-time function, x(t), is sampled every DT time 
units, the resulting samples have values x(n-DT). These samples 
produce a sequence of values denoted x(n) as used in the preceeding 
chapters, The values can also be denoted x and used 7 define a 


vector x: 
EC o O M 


The sequence x has a Discrete Fourier Transform, denoted x: 


N-1 E j2mkn 
- — 46 N 
^ n 
n=0 


The values of the DFT can also be represented as defining a vector 


ам 


x= (xxx yA 
25 202512“ | 


ln matrix notatation, the DFT relation can be expressed as: 
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UR x 
Or; 
1 
Xo 1 1 557”. Xo 
2 N-1 
X4 1 W W “ > 7 PW X, 
2 4 2(N-1) 
(5,1) xə d 1 W W ə cc W MEL 
: : Q1) 20-1) y (1-1)? 
"IET LE LAN 


where W= exp(- ST), 

It is desired to trace the effects of the noise from the time 
domain to the quefrency domain, The first step therefore is to 
examine the behavior of the DFT of a noise Sequence, To facilitate 
the derivation, the noise sequence to be encountered will be considered 
to consist of samples of a zero mean, complex, white gaussian random 
process, Ihe appropriateness of this model is 5000507 in Appendix D, 


To proceed, make the folloving definition, 


[> 
ii 
I< 


v is a zero mean, complex, white 


gaussian random vector 


4. 


The sequence x then has the following properties 


E(x] = E[r] = £ İlm(ə İ - o 





-80- 


E [re Re(x)| = E: Is I = Nx N identity matrix 
в | ты?) n(x) | = E I 

E | ne( ) таб) | = О 

Е İz” xs) uz T where * denotes complex conjugate 


In short, x is a vector of N complex random variables. The vector 
contains a total of 2N uncorrelated zero mean gaussian random 
variables (GRV's). The uncorrelated property, for zero mean GRV's 
implies independence. The complex variables, xə can also be represented 
as, 

HS Ix, |-exp [ARG(x, )] 

where Ix, | is & Rayleigh distributed random variable and ARG(x, ) is 
a random variable uniformly distributed on (-4,17). 

One way to arrive at the statistics of x is to argue along 


the following lines. From equation (5.1), 


2k ‚k(n-1) 
+ ut X4 + W xə s T y Xy 


4 
| 


Ak XQ 1 


| 


9 ¢ 9 § 
(5.2) xo + X4 + xə + cee t+ ХЫ 4 


иһеге x] = Wee s The purpose of defining this new sequence x. is 


as follows: 
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Ix: = ГБ x, | = we] 7. = |x, | 


|x; [is therefore Rayleigh distributed with statistics identical to 


those of |x; , Furthermore, ARG(x; ) E ARG(V”T) * ARG(x, ). Since 





the phases are computed as principal values, ARG(x;) Toma So) unl 
formly distributed on (1,1). Consequently, X: has Rayleigh magnitude 
and uniform phase values which are SUE identical to those 

of X; It follows that the real and imaginary parts of x; must also 
be statistically identical to those of Xs i.e., independent zero 


mean, identically distributed GRV's. It is also true that the 


elements of the x' vector are independent: 


E E xt = İn [wm ә x 
nen : m 


| 
m 
in 
s 
1 
з 
=” 
zi 
Ы 
s 
ə 
* 
SS Ј 
l 
a 
ә 
О” 


By the zero mean GRV property, this implies independence, The result 


of this intermediate step can be seen by re-writing equation (5.2) 


as below. 
(5.2) Xx = Xo T хј zi xə Ti. T Xx 

Re(x, ) = Re(x$) + Re(xj) 057 Re(x& 4) 
(5.3) 


In(x,) Im(x$) + In(x}) Boa... + In(xy 4) 





E 
The right sides of equation (5.3) contain a total of 2N independent 
zero mean GRV's. The two sums, i.e., R(x, ) and In(z,), are therefore 


independent. Additionally, they each are sums of N independent zero 
2 + 
— Hü Tolloysahen that Re(x,) and 


In(x,) are independent zero mean, identically distributed, GRV's with 


mean GRV's with variances 


2 
ES 
2 


It is also possible to show that the individual elements of 


variances N so that is Rayleigh and ARG(x,) is uniform. 
Y, k 


the x-vector are independent. This follows from the orthogonality 
phe rowsoot the Wzmatrix, 1.e., Др and x, become orthogonal 
random variables. Since they are also zero mean and gaussian, they 


are independent. To see this more explicitly, consider: 


E ls, =) = E b .—5.x. eg 


A -m -m(N-1) 
kz + W Xİ us M 4а 


Since E x x | = o” 05:04: Tollovslhat 
157: ij 


R |l x y(k-n) bx ien Den). P 


tz 
ә 
V 
Fick 
C 
il 


e — yn n) 


n=0 


For k = m, the value of the above sum is clearly No”, For k = m, the 
sum consists of N terms, all with equal magnitudes and with phases 


equally spaced from -m to Te The resultant sum is therefore zero so 


that, 





N-1 .: 


E Ек m) = ос ¥ unm) AW ö km 
n=0 


Once again, together with the zero mean gaussian property, this implies 
independence. 

In summary, a zero mean, Complex, white gaussian noise sequence 
has a DFT which is also a zero mean, complex, white gaussian noise 
sequence, All that changes, statistically, is that the variances 
are changed from 2 to No”, 

In computing a DFT, it is sometimes desired to add zeroes to the 
end of the sequence, Since the same DT is used, this does nothing to 
change the bandwidth of the DFT. However, since there are more samples 
of the time sequence, the DFT will have more samples in it. This has 
the effect of increasing the resolution of the DFT. The effect of this 
operation on the statistics of the noise DFT should be examined. 

suppose that M samples of the noise sequence, x» are available 


and that an N point DFT is computed: 


= 1 i E c M дә 
= 1 Y we mər x, 
xə 1 We W yə) X, 
Да : ə Ф ° 

. > : ə 
ı “a. : 0 
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The x-vector will, in general, have N non-zero elements which, 
as above, can be shown to have real and imaginary parts which are 
independent, zero mean GRV's with variances M :7 Hence the 
magnitudes of each point are Rayleigh distributed and ihe phases 
are uniformly distributed. The sequence, however, is no longer 
white. In fact, the x-vector Will span the same M dimensional 
complex (or 2M dimensional real) space as that generated by the 
original sequence x. The autocorrelation function of the DFT can 
be computed in a straightforward manner. From this it can be 
shown that, for example, if N/M = 4, points in the DFT 4 samples 
apart are independent. If N/M is not an integer, all the points 
are correlated, although only loosely for seperations of more than 
about 2N/M. 

If X =s + v, where s is a deterministic vector and v is the 


noise vector as previously described, it follows from the linearity 


of the DFT operation that 


AS E [Sk + xk] = Bk 
and 
i ES xx] - E (A E NC. | 
_ 2 
~ "Se Sn Oe 


With this background, consider the DFT of the sequence 


Xa = S, + V, Which may appear as in figure 5-1. 
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Figure 5-1 Frequency Plot of Signal and Noise 


The solid line in the figure represents the magnitude of the 
DFT of the deterministic sequence. The dashed line represents some 
convenient probabalistie bound on the magnitude of the DFT of the 
noise as determined by its Rayleigh distribution. What must be con- 
sidered is the effect of tne noise on the phase of the composite 
signal. For values of k such that Is, | > Iv |ә the vector sum may 


be as shown in figure 5-2(a). 





Figure 5-2 Vector Diagram of Signal and Noise 
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In this case, the composite phase will follow that of s, reason- 


k 
ably closely. As an example, if E = 2 |Va. |» the phase error, the 
difference between are (s, ) and are (s, T Y.) is at most about e TM The 
phase curve is therefore noisy; but the noise is small - less than 

1/2 in magnitude. 

For values of k such that İkİ « Iv.) the vector sum may be as 
shown in figure 5-2(b). The phase error in this case can be considered 
to produce à more noisy curve - errors as great as m in magnitude. 
However, the implications of this situation are more significant in 
view of the fact that the phase curve will now have the structure of 


the noise phase, i.e. white. This causes problems because it is 


presumed in the phase unwrapping algorithm that the phase has a smooth, 


or continuous, structure. If BA > [y |for all k, every point in the 
DFT phase will be within +5 O2 theseorreet value so that the measured 


number of branchings of the arctangent function defining the composite 

phase will be, at most, one different from the correct value, If 

Ix | « Wich however, there can be any number of branchings of the phase 

So that a huge error can be expected in the corrected phase curve. 
This can be summarized by considering the form of the complex 


log of the DFT. 


log e + Y) = log Әһ * log 11 7 
Y 
s logs, + — ? lex | >> xx! 


vk 


LX 
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S 
A: 
log PESA = log Yg t Tog m 
Pk 
2 — , > 
~ log Y. + Y. MN > EN 


y 

If the phase of y, is uniform and white, the phase of = will also 
k 

be uniform and white, when measured via its principal value. Therefore, 
Yg Re(y,) + j In(y,) 
28 statistically equivalent to ua Du Consequently, 
^k mk ; 
the effect of the log DFT operation for ЕР” [ys] is Loss rcr tne 


noise as in figure 5-3(a). The |s,|used is that of the Loran-C signal. 








(b) ? 


Figure 5-3 LOG DFT Filtering Effect B 
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For BA < [Yk |» the composite log DFT is essentially that of Vie 


Consequently, the resultant log DFT can be approximated as: 


log (Sy * X = [108 sk) OH 4 108 pee i —. 


where H, and H, are as shown in figure 5-2 (Dana (с) and x | and 


2 
ARGÜYL) are the Rayleigh magnitude and uniform phase as previously 
derived. The equivalent operation is depicted in figure 5-4, 

The implications of this problem can now be seen. The output of 
the phase unwrapping block of figure 5-4 will be grossly in error if 
the Ko of figure 5-3 corresponds to a frequency less than the sampling 
frequency, 1/2DT. For a given signal to noise level therefore, sampling 
faster than a prescribed rate has dire consequences. It is, however, 
desired to sample the input signal as fast as possible, obtaining as 
many samples of the groundwave as possible before the first skywave 
arrives. The number of samples so obtained will define n, in İhe 
complex cepstrum, the quefrency seperation between the BLA terms 
of $ (n). If this number is reasonable, e.g., ^ or more, the estimators 
can be expected to perform well. If it is less than 4, there will 
be too much overlap between the terms of P. (n) to obtain the desired 
precision, If a greater bandwidth is necessary, some sort of 
smoothing would have to be applied to the log magnitude and phase 
curves to reduce the errors introduced by the false branchings due 


to noise, In such an operation, however, precision is lost in the 


cepstrum, 
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Consequently, a search for a compromise along the following lines 
must be made. 

1. For high signal to noise ratio, determine the ranges 
of values of sampling rate and arrival times in which 
acceptable performance of the cepstrum estimators is 
achieved. 

2. Corresponding to values within these ranges of sampling 
rate and arrival times, determine the lowest signal 
to noise ratio for which the estimators will function. 

The minimum SNR will determine how many pulses will have to be 
sampled to obtain the desired performance. Through considerations 
such as signal stability, a means of evaluating the technique is 
available. 

An example should now be considered to see how well the theory 
developed thus far fits the results. Signal power is typically 
defined in terms of the power at the 25 microsecond sampling point. 
As it turns out, this would be consistent with most other definitions 
of average power, Consider sampling every 10 microseconds for a 
500 microsecond duration and then taking a 512 point DFT (N/M = 10). 
Here, the normalized signal power is, s_ ? « 0.26. For each 
component of O db noise, the variance is therefore, 5 = 0,13, For 


+40 db SNR, g^ - .000013. In the DFT, y 


- 50 g^ = 0.00065 so that 
g = “0245, 

From the Rayleigh distribution, it can be determined that 
Sk | > Мә with confidence levels „99 and „86 for flog Er x -3.0 


and -2.6 respectively with v as above, From examination of a noise 


k 





20 
free DFT it can be verified that these values are achieved at 
frequencies of about 27 kHz and 31 kHz. The case was simulated and 
the resulting magnitude and phase curves are presented in figure 5-5. 
From the phase curve it can be seen that the noise induced branching 
begins alone 36 kHz. 

For a +30 db SNR, the same confidence levels correspond to 
frequencies of about 17 kHz and 19.5 kHz. The results of a simu- 
lation are shown in figure 5-6 where it can be seen that the first 
branching occurs at about 23 kHz. For +20 db, the predicted values 
are 11 kHz and 13 kHz while, as shown in figure 5-7, the branching 
occurs at about 18 kHz. If the magnitudes of these plots are 
considered, it can be seen that the frequencies at which branchings 
could have occurred are nearly exactly as predicted. 

The cepstrum for the +40 db case was computed with the following 
Signal parameters: f. = 0, DT = 20 psec, To = 20 usec, Ti -To = 80 usec, 
B, = (3,3). The resulting cepstrum is shown in figure 5-8, As can be 
seen, the error is small so that the estimators can work properly. As 
an indication of this, the estimate of the groundwave phase was 
-0,00822 radians, about a 13 nanosecond error. 

In the +30 db case, the phase unwrapping presents problems since 
the 25 kHz point, i.e., the k = N/A point, is above the frequency 


at which İYİ NI e The cepstrum is shown in figure 5-9. 





These examples are felt to verify the theoretical limitations 
of the technique as presented earlier in this chapter. A discussion 


of the implications is included in the following chapter. 
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Figure 5-5 + 40 db SNR LOG DFT 
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Figure 5-6 + 30 db SNR LOG DFT 
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Figure 5-7 + 20 db SNR LOG DFT 
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Figure 5-8 + 40 db SNR Cepstrum 
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Figure 5-9 + 30 db SNR Cepstrun 
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Chapter 6 


CONCLUSIONS 


The results of the preceeding chapters will now be summarized, 

The cepstrum of an undistorted received Loran-C signal consists 
of components which are known to within & few parameters. If the 
noise level is such that & sufficient sampling rate can be employed, 
the following condition ensues. While there is overlap and interference 
Berveenzise componentsgeinsthiescepstrun sthererare rangcssetequerfrenesress 
albeit small, in which particular components dominate the sum. In these 
regions, parameters associated with the dominant component can be 
estimated. 

The performance of the estimators under ideal noise conditions 
was demonstrated in Chapter 3. The technique estimates the groundwave 
phase as well as could be desired and has the added advantage of 
producing estimates of the various other parameters of the signal 
which could be of interest. 

The results of Chapter 4 show that the signal distortion intro- 
duced by the propagation paths, of the amount to be anticipated in 
the Loran-C case, has negligible effect on the estimators. It was also 
shown that the distortion produced by the receiver filters could be 
accounted for so that performance is maintained - if the filter 
responses are reasonably well known, 

The results of Chapter 5, although not conclusive, give an 
indication of the amount of processing time required by the scheme, 


In the +40 db SNR case, a bandwidth of about 36 kHz was usable. A 
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corresponding sampling rate of about 14 microseconds would produce 
values of about 3 and -.5 for ny and A, respectively in the case of 
a 35 microsecond skywave arrival time. Such a value of n4 is a 
limit for the type of estimators used here. This is due to the fact 
that the A, and B, estimators can not make measurements on y(n,-2) 
if ny = 2 since the zero quefrency point is where all of the 
groundwave magnitude and phase information is mapped to. 

It has been the purpose of this thesis to consider only the 
computational aspects of the application. Consequently, a conclusive 
statement 75: , of implementing the technique can not 
be made, What can be done however, is to cite some practical 
limitations of the Loran-C system and to review the results in 
ligni of inem. 

The major conclusion of the results is that a signal to noise 
ratio of about + 40 db is required. If each individual pulse was 
at a zero db power level, 10,000 pulses would have to be averaged to 
achieve the desired performance level, If all eight pulses of a 
group with a 100,000 microsecond repitition interval were used, 
this would require an averaging time of just over 2 minutes. Such 
long times require great signal stability and, of course, preclude 
use in a navigational application, i.e., on an aircraft or ship. 

The stability could be achieved in a stationary application, as in 
a calibration or monitoring receiver, but other effects must be 
considered, 

The eight transmitted pulses are not actually identical. 


There 1s some degree of droop in signal power since it is difficult 





mI. 
for signal transmitters to recover from the effects of a large power 
pulse generated only 1000 microseconds earlier, The effects of 
such droop can be taken into account somewhat in the sampling- 
detection schemes presently employed. For a homomorphic filter 
receiver,however, it may be required that the averaging be accomplished 
on the first pulse of each group, thereby making the averaging time 
16 minutes, 

In addition, there are the effects of non-linear filters used 
to suppress the non-gaussian aspects of the atmospheric noise. These 
and other non-linearities involved in actual processing would have to 
be taken into account, especially with regard to whether they change 
the nature of the signal drastically. 

Speculation on the feasibility of the technique therefore leads 
to the following conclusions. The technique shows good results 
when applied to the most general signals that can be simulated. In 
view of the long averaging times and the more complex nature of an 
actual implementation however, the improvement in performance does 
not seem marked enough to conclude that it could be effectively 


implemented. 
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Appendix A 


COMPLEX CEPSTRUM DERIVATIONS 


This Appendix presents derivations of the complex cepstra of the 
various signals considered in this thesis. The complex cepstrum of 
the Loran-C signal, in its most general form, is somewhat complicated, 


It can, however, be arrived at in the following systematic fashion, 


A-1 Basic Waveform, Arrival Time Coincident with Sampling Time 


Consider first the basic waveform, s(t) = 5. и (0), sampled 
so that: s(n) = ne le? u qn), The first step is to compute the 
Fourier Transform of this sequence. It can be most easily calculated 


by utilization of the following property. 
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then, S- P(e” 
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and, 
dw n =-00 


The term on the right side of the above equation is, by definition, 


the Fourier Transform of -n^f(n). It follows then that, 


2 | a“ jw 
n”f(n) $— - =o She ) 
aw 


To utilize this property, let f(n) = e" u_,(n). Then, 


iw —.. -јип -(ал-ји) E 
F(e!") = ie ју". zc DD WADE. 
n=0 
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and, 2 s 2 E : _ 
Eey o -fa loo] 
dw dw 


on (arin) fy, late] 
we -= O4-jw | E 





which is, therefore, the Fourier Transform of s(n). 


As described in Appendix B, the computation method used in the 


implementation will remove the linear phase tern, ence thereby yielding 


a computed transforn, ( 
: 5 - (Atjw 
(ed) " a ite 


The next step is to compute the complex natural logarithm of 


ste, 











ј - (A+jw T 
log S (eJ") = -Q + Log] 14 C (ол | 516 I Tr) 
Since A >0, and therefore, (rm) < 1, the following logarithm 
expansions are appropriate: 
n T | 
logli+x) = E QUIE, lx|< 1 
n=1 > 
co x 
106(—-—) 2 n : ҝ|- : 
-X 
n=1 
Therefore, 
co m. -0n kİ 
log S(e!") = -Q+ 2 co" ER c јын) 
n=1 
(A.1) È 
eS -an : 
+3 E | 5 
n=1L 2 


The final step is to evaluate the inverse Fourier Transform of 


log S(eJ") anà to denote it as £(n), the complex cepstrum of s(n). In 
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this case, the form of the expression in equation (A.1) allows the 
inverse transform to be obtained by inspection since (n) is the 
JWN 


coefficient of e in its transform. 


(A.2) S(n) = -a — 


nti 
get (=) ae end ng 


Two generalizations can be made on the results obtained thus far, 
For any constant, B, 


jw 





B s(n) B S(e"'")——— log B + log ee) orqa B)ö(n) + S(n) 


so that any coefficient change, real or complex, in the original signal 
effects only the n — 0 point in the cepstrum - in an additive fashion. 
Additionally, consider the effects of sampling the continuous-time 


function, s(t) = ак" u 4(t), at a rate (1/DT). This yields, 


(n: DT)? ,"O(n-DT) 


s(n) 


İl 


u_, (n-DT) 


(DT)? „2 a "NO 
Thus, there is a coefficient change and a virtual change ing. 
Therefore, to compute the cepstrum of s(n) for various values of DT, 
simply modify the expression in equation (A.2) by changing the exponent 


from a to (a-DT) and adding the term (log DT) 6(n). 


A-2 Basic Waveform Arrival Time Between Sampling Times 


The notation in this case becomes a bit awkward and requires explana- 
tion. The expression ö(n-x), for non-integer x will be used. What is 


implied is that s(n) GQ 6(n-x) is a sampled version of s(t) x 6(t-x). 
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With this notation, consider 
s'(n) - s(n) ® ein) 
where g(n) = 6(n-A) 2001777770 
S'(eJ") » s(eJ") c(eJ") 
S(e7") will be as derived in the preceeding section and 
G(eJV) = 9 UWA 
log G(eJ") = -3wA 


The contribution to the cepstrum of s'(n) due to g(n) is, 


. T : 
@(n) = p [105 ofe") | = = -jW À ev” ay 
ili 
— aon jun (j T - jTn ә 1 | 
= z ļe jm - 1) +e (јал = 1) 
2Tn = 
— uu... HR. 
-— 2 23 


2 2 


A | sin nt m cos 4 
1T 
n n 


For n - O, 
Дур ә TL y Ди 
e(n) = = fo əl on |= (-1) = 
For n = 0, consider n to be a continuous variable and consider the 
limit, 
lim ES m__- m cos = 
2 
n— 0 n 
13 = cos nil = T cos ПП 7: sin 4 
= inm — ——— —— — - ——-. —.-. ———-.—o— 
0 2n 


lim E sin s = 0 
иә 
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Consequently, 


- (aye A azo 


Since s'(n) = s(n) 8 g(n), so that $'(n) » S(n) *'/E(n), it follows 
that, 


) n+1 EA 


(-1 à 


25250 
S Ww = -Q none) 


n+1 n+1 
CFE gt Gl an XP 
n n 


A-3 Multipath Impulse Train 
The impulse train, denoted p(n), consists of 
(2:3) p(n) = 6(n) + B, 6(n-n, ) + t B, ӧ(п-п ) 


To begin the derivation, consider the case of m = 1, [B4] =], 


p(n) = 6(n) + B, 6(n-n, ) 


oO 


P(ed") = = p(n) e dn 5 edəni 
==00 
log P(eJ) = log | 1 + B, eJ | 
E k 
® (B, en 
= > ae. 1 since İB < 1 
ket 3 : 


The cepstrum is then, 


Sn) = FT E prod) | 


n 
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k 
co B AT -ju(n-kn, ) 
E) = E» = il —һӧм- ди 
К=1 -Ti 
co k : Tf 
" Y (yt A í,"Jw(n-kn, 
k-1 k j2n(n-kn,) -T 
k 
oo B sin t(n-kn, ) 
N 3 = E: k+1 m 1 
(A.4) Din) " (-1) k “vü 


If n, is an integer, equation (A.4) becomes 


k 


co B 
Sn) = zr [C07 -i 6(n-kn,) 
kei 


so that the echo contributes an impulse train to the cepstrum (for 
integer пу). ln general, of course, the more involved ird behavior 
of equation (A.4) is encountered. 


The more general form of equation (A.3) calls for more than one 


echo. In the case of m-2, 


p(n) = ӧ(п) = B, 6(n-n,) * B, ö(n-n,). 
In this case, 
P(e") = 1+ Beyl + BD 


If EE İB. € 1, then 


log P(e") log E per se + xu | 


Zu = 
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(B e? 1 4 B.e 
en il 7 2 
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This is the inverse Fourier Transform of the cepstrun, 


2 2 
: . B B : 
ЈА, АТ О - jim à 1  -jw2n ieee JWOD 
F9] = B.e 1 > E 1 + ə e 1 = eee 
-jwn P -jw2n DT - ju3n r 
+ B.e 2 - -5 e 2 + ns e 2 
“nuz EHE En + B, e у"т?) 
ve 
. p,p,20"31(n,+2n,) 2 par An oe 


By a straightforward extension of the previous results, it can be seen 


that the cepstrum will consist of samples of -. functions 


2 2! e... Кпдә 


functions centered at all 


centered at quefrencies nas en vo kn, and n,, 2n 


Tu 


Additionally, there are samples of —— 





combinations of multiples of n4 and nəə 
Clearly, as the number of echoes increases, the cepstrum of the 
impulse train quickly becomes complicated and difficult to filter, To 
overcome some of the complications however, use can be made of the 
form of the signal information transfer between the time and 
quefrency domains, As Sr has shown, if the cepstrum of a 
signal is short pass filtered with a cutoff quefrency of, for example, 
Mos the recovered waveform is identical to the original waveform for 


time values up to m Additionally, for reasonably large values of m 


0 0” 
the recovered waveform follows the original signal very closely for some 
time after mo" Due to the pulse nature of the Loran signal, this fact 


makes a combination filter effective, 





-108- 

The cepstrum can be short pass filtered with a cutoff quefrency 
of about 150 usec, This will retain almost all of the groundwave 
information butremove much of the complicated impulse train contri- 
bution, A simple example can illustrate the effectiveness of this 
approach. Suppose n, s ДО usec, nə s 80 usec, and nəs 120 usec. 
If a cutoff quefrency of 150 usec is used, sampled SIRE functions 
centered at the following quefrencies e reiained: 

D, 2n, 2n, , un, 


nəs “nə 


= 

n4*n,, “ni nəə An m, 
Some of the above quefrencies are above the 150 usec cutoff, but close 
enough that their effects may be felt in the pass band. Other sampled 
sins functions are centered at quefrencies far enough removed from 
the cutoff that they may be considered removed by the short pass 


filter. Consequently only the quefrencies listed above have to be 


considered in the more elaborate filter-estimators of Chapter 3. 


A-4  Cepstrum for Exponentially Weighted Signals 


The multipath impulse train thus far has been considered to 
be 
p(n) = ö (n-no) 25 6 (n-n,) + B, 6 (n-n,) teos +B ö (n-n,) 


with the following condition met: 


k 
(4.5) — i İb E 1 
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This ensured that the groundwave was the dominant component of the 
composite signal and allowed the expansion of log в(е9") in the 
manner presented, As will be shown in Appendix B, this is part of a set 
of conditions sufficient to ensure a complex cepstrum of the form 
desired for the methods of this thesis, 

In the interesting Loran-C cases this condition is typically not 
met. The problems can be circumvented however, by exponentially 


weighting the signal in the manner to be described. Suppose 


x(t) 


s(t) + B, s(t-T,) 


mE cct 2 -a(t-T, ) 
= te u,(t) + B, (t-7,)° e 1321) 


vith EA = 
| -Bt : 
If the signal is multiplied by e” , vhat results is 


Qt x(t) = te (08)t u(t) + B (t-T, )* o7 (048) (t-14) u_, (t-7,) 


where B, = B, er 
If B is such that [B; | < 1, the condition of equation (A.5) is 

satisfied, The resulting cepstrum will be as previously derived with 

a replaced by (at8) and B, replaced by Bj. The cepstrum can then be 

computed and filtered in the manner developed in this thesis, If it 

is desired to recover the original signal after filtering, the filtered 

+t 


signal can be multiplied by e as a final step. 
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The results of the above argument can be generalized to the 


case of multiple echoes to obtain the requirement that 


k 
» |в je" "m < [Boje 7o 
m=1 
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Appendix B 
PHASE UNWRAPPING CONSIDERATIONS 


As mentioned in Chapter 1, the computation of the complex 
Dearth of the DFT of the signal presents problems. A computation on 
a point by point basis via an arctangent power series expansion will only 
provide the principle value of the desired phase, What would result 
in the case of a typical loran-C signal might produce a phase plot as 


in figure B-1 (the dotted lines). 


- —:— ad 
. » N/2 
=T 
-2m 
-Зп 
T: 


Figure B-1  Uncorrected Phase Curve 
It is desired to remove the discontinuities induced by the branch- 
ing of the arctangent function for the following reason. If x(n) = 


s(n) © p(n), then it is desired that 


log x (e)? - S(e?") + log P(e") 
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In particular, the above equality must hold for both the real and 


imaginary parts of the expression, Thus, 


log X(eJ") = log S(e"”) + log P(e”) 
(3.1) and, 
atg X(e7") - arg S(eJ") + arg P(e7") 


If, for example, there is some value of w for which arg S(eJ") 
has a value 3n/4, and arg P(eJ") has a value 7/2, then the right side 
of equation (B.1) equals 54/24, However, arg He), as evaluated by 
the arctangent power series will have a value of -31m/4 so that the 
equality is violated. 

SX has described a method for removing the discontinuities, 
After the complex log of E DFT of the signal is computed, a sequence 
with imaginary part appearing as the dotted lines in figure (B-1) 
results, The imaginary part of the sequence is then scanned to detect 
jumps of almost 2T in consecutive values, For appropriate values of 
k at which these jumps occur, a correction sequence is generated. 

This is represented by the heavy lines in the figure, The correction 
sequence is added to the originally computed phase sequence as shown 
in figure B-2. 

As = be seen from the corrected phase, there is typically a 
large linear phase term involved, Since this term can easily dominate 
the phase plot, and the resulting complex cepstrum, it should be 
removed if it has no effect on the results, Schafer's method is to 


examine the value of the corrected phase curve for values of k near N/2. 





lo 





Figure B-2 Corrected Phase Curve 


The integer multiple of q that is closest to this value is called the 
linear phase term. From this value, the straight line shown in figure 
B-2 can be generated. If this is subtracted from the corrected phase 


curve, a final version of the phase, as in figure B-3 results. 





Figure B-3 Corrected Phase Curve - Final Version 
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The effect of this last operation is to multiply the DFT of the 
signal by ERU where mə is the integer linear phase term deduced in 
the manner described above, This is equivalent to convolving the signal 
with 6 (n-m,) - or shifting the sequence mo time samples, As a final 
step after the filtering is accomplished and the inverse operation 
is performed, the recovered signal can be re-shifted this amount to 
produce an output with a proper time base. 

The corrected complex logarithm is now continuous and periodic, 
and the desired complex cepstrum can be computed from it. In effect, 
what is presumed here is that the branchings of the phase curve are 
caused, in a systematic fashion throughout the complex log sequence, 
by the same linear phase term. This presumption should be examined 
more closely. Consider the impulse train component of a typical 


echoed signal, 


k 
mn. nn) 
i-o 3 : 


The Fourier Transform of this sequence is 


i k , . 
PME B. e Joni where B, - İB, (67: 
1-0 


The following analysis will consider the form of Р(2У") LOr 
various values of the İb, İs, 2,*5, and n, "s, To begin the analysis, 


= 0, andk -1, Furthermore, let B, be real, 


let By = 1, 2, = О, No 1 


positive and less than 1 with n, = an integer, What results in this 


case is, 
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P(e2") = Pl o By su 
n, = an integer 


A vector diagram representation of P(e") is as shown in figure B-4, 





Figure B-4 Vector Representation of P(e7") 


Às w ranges from -m to mq, the term DK traces out the circle 


n, times, If B, is complex, i.e. B, - |B, | oy, 9, 2 0, the dashed 


line of figure B-4 begins (for w - 0) at an angle 9, with the line 





representing the vector 1. Importantly, as long as By sm ihe branch 
point at «+1 of the arctangent function defining the phase is avoided. 
Consequently the magnitude and phase of P(e3") are as plotted in figure 
B-5, | 

For |, | « 1, A, the maximum value of the phase in figure B-5(b), 


is always less than 71/2. ar m is not an integer, the Be 974 vector 
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ARG P(eJ") 


NAN 
AINE V 


Figure B-5 Hagnitude and Phase of p(e%) 





in figure B-4 starts in line with the 1 vector but ends up, for w=", 
at some arbitrary angle. In terms of figure B-5, this means the plot 
is essentially the same - except that the w = HT re are no longer 
nulls or peaks of the magnitude curve, nor zeroes of the phase curve, 
For B = İb, İ e, 2, 2 0, the entire plot of figure B-5 is shifted 
to ihe right or left, 

With these preliminaries, more general impulse trains can be 


considered, Suppose 
p(n) = ö(n-n,) T By ö(n-n,-n,) o 


Consider first, the case wherein ng 7 I; B4 = a real, positive 





SI 
nunber less than 1. The magnítude of Plet") will be as in figure 
B-5(a). The phase however, will be as in figure B-6(a). When the 
two step correction procedure is employed as previously described, what 
results is presented in figure B-6(c) and is seen to be identical to the 


phase curve of figure B-5(b). 





Figure B-6 Frequency Representation of P(e") With Linear Phase 
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A vector diagram of P(e”), garametricoln v, isfshoövnutor this 
case in figure B-7. Notice that the total resultant number of crossings 


of the branch line is one, i.e. n This is true independent of the 


“oie 
orientation of the axes in the figure, i.e. for arbitrary By 


= 


e [»(e7*] 


= 


m. = i 


Figure B-7 Vector Representation of P(e”) with Linear Phase 





The following summary can be made of the results thus far. If 


p, (n) m ӧ(п) * B4 6(n-n, ) |B; | <1 





— 


and, 
pə(n) = ö(n-ng) t B, 6(n-n, -ng İB < 1 
e 
integer 
then, pə(n) = p, (n) x ó(n-ng), 


imo. pj (n) is p, (n) shifted ny samples to the right. Since any linear 
phase term is removed in the computation, p, (n) and p^ (n) become 
identical signals beyond the complex log stage in the processing, It 
need merely be remembered to shift the recovered signal back no 

Samples as a final step in the second case. 


Now consider the case wherein 


p(n) = ö(n) + B, ö(n-n,) By = |By| > 1. 
The resulting phase curve is as shown in figure B-8. What can be 
seen is that the linear phase term is due to D. The ripples in 
B-8(c) are as in figure B-5(b) but reversed in sense. This can be 
 iiributei to the fact that the 1 vector in p(eJ") is rotating 
counterclockwise relative to the linear phase term associated with nye 
Further light can be shown on the effect by considering the logarithm 


expansion of P(e") in this case. 


log P(e%”) 


İl 


log E + By ml 


- loc o, emi 577. |: T в 7 — 


-k 

co B : 
= -jm * logB, + 2 ei, E edwin; 
k=1 





(a) -120- 


(b) 





Figure B-8 P(e") Phase, B4 = 1 


— 
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With the linear phase term removed, the cepstrum becomes 


-k 


oo B aT : 
$() - Qog B) 6() *  z(-D M Hs mini) ay 
Кеј E 


n 


co kt+i B 
(log B,) 6(n) + BE — банкы 


for an integer Nye This means that the impulse train appears in the 
negative quefrency portion of the cepstrum. What is further implied 
is that the largest component of the signal contributes to the low 
quefrency portion of the cepstrum while the smaller components 
produce the impulse trains. 
In the Loran-C case, it is desired to recover the groundwave 

(or at least to treat it as the primary component). Consequently 
exponential weighting is used to ensure that the groundwave is the 
largest signal. This also simplifies the filtering process since, in 
a multiple echo case, in which a middle echo is the largest, Be ee 
would appear in both the negative and positive quefrency regions, 
Additionally, various cross product terms in the logarithm expansion 
Would produce impulses throughout the cepstrum (e.g., at n, - 2n,, 
nz - Jn, , 567 

From figure B.8(b), it can be seen how easily the linear phase 
component can obscure the important information in the phase sequence. 
Clearly, it is important that this term be removed. SE, in a 


multiple echo case, the discontinuities need not be caused merely 


by the linear phase component of the dominant signal. In terms of 





PV 
the logarithm expansion, 
jwn 


: k 
log Ple") = log | X B. dE i 
i=0 


the requirement can be seen to be 


- wn. 
e J i 





k 
2 B, 





(B.2) [Bol > 


to ensure that branching is due only to the linear phase term 
associated with E 
Schafer states that, after low pass filtering to avoid aliasing, 


the condition, 


k 
(2.3) 2 


is sufficient to meet the requirement of equation (B.2). This is 
certainly consistent with everything presented thus far and would be 
true if an ideal analog low pass filter were employed. In such a 
case, the condition of equation (B.3) could be met via exponential 
weighting as presented in Appendix A. However, for an analog filter 
with a finite passband to stopband transition width, the effects of 
aliasing should be examined more closely. 

As it develops, for transforms with ripples as in figure B-5(a), 
the effects of aliasing can cause large amounts of distortion in the 
phase which can not be overcome by analog filtering. To show this, 
a more in-depth development of the DFT must be presented. As an example, 
the DFT of a sequence is represented by the dotted lines in figure 


B~9 (the sequence is e u 4(n) ). The solid lines in the figure 





Ed 


kedal Part 





Imaginary Part (b) 


Figure B-9 Real and Imaginary Part of le DFT 


represent the continuous Fourier Transform of the corresponding 
continuous time function depicted in a manner appropriate for compari- 
son with the DFT, i.e., with the negative frequencies shown to the 
left ofk = N. 

It would, of course, be desired that the DFT would consist of 
samples of the appropriate solid lines, i.e., samples of curve 1 for 
k « N/2 and samples of curve 2 for k > N/2. As is the case however, 
the DFT is actually samples of the sum of the two curves for all 
values of k, Due to the bandwidth of the signal and the sampling 
rate employed in this case, it appears that the aliasing has a negligible 


effect. The only errors that appear remotely discernible are in the 
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imaginary part of the DFT for values of k near N/2, and even these 
are small. 
If the DFT is shown via its magnitude and phase as in figure 
B-10 however, the evaluation changes. There is a large error in 


the phase for a wide range of values of k in the vicinity of N/2. 






Magnitude 





Figure B-10 Magnitude and Phase of e u 40n) DFT 


What can be seen is that for k < N/2, the composite phase is 
in the same quadrant as that of curve 1. For k > N/2, the composite 
phase is in the same quadrant as the phase of curve pu This is due 
to the fact that the magnitude of the continuous Fourier Transform of 
-t 


e u_,(t) decreases monotonically with frequency. Consequently the 


aliased signal (curve 1 for k > N/2, curve 2 for k « N/2) is always 
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smaller than the desired signal so that the sum is dominated by the 
correct signal and the resulting phase is of the correct sign. 

With this background, the actual signals treated in this thesis 
should be considered to show what can possibly go wrong. Suppose 
x(n) = s(n) @ p(n) 
so that 


X(eJ") - s(oJ") p(e3") 


Iset]: fec) exp | Jj E ә(е9") + arg ple") 





Magnitude 


yed 


Phase 
31/2 


(——————— - 


-Зп/2 1. | 


Figure B-11 Magnitude and Phase of Loran-C Signal Transform 





5757 
The magnitude and phase of the continuous Fourier Transform of s(t), 
the Loran-C signal, are shown in figure B-11. What is important 
to notice is that, for values of k less than, but nearly equal to, 
N/2, curves 1 and 2 are nearly 180° out of phase (modulo 21), 


Now suppose that 


p(n) 6(n) + B, 9 (n-n,-A,) 


i 


where B, 0.8 
ny = an integer 
A, Es -.5 
so that 
P(e") = 1 4879" (пун) 


mesvector plot of this Р(е9") is as shown in figure B-İZ, 


ӧн | 1.8 
ine 
line 3 = 
и = (п-ӧ) 


line 5 

w= -(7+5) 
line 4 
W = -T 





Figure B-12 Vector Representation of P(e?) 
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The rotating vector begins, for w = O0, at the position shown as line 1, 
As w increases to tm, the rotating vector moves in a clockwise 
direction to the position shown as line 2, When w is at a value of 
+(1-5), the vector is at the position of line 3 - where the resultant 
magnitude of р(е9Ҹ) is 0.2. As w decreases to a value of -T, the 
rotating vector moves in a counterclockwise direction to the position 
of line 4. As w continues to decrease, the vector arrives at the 
position of line 5 when w = -(m+9) - where the resultant magnitude of 
P(e") is 1.8. 


The implications of all of this is as follows, let 


X Se KX FM where X = DET 
^ 1 2 дә 
X, = desired 
1 : 
signal 
X. = aliased 
2 : 
signal 


Then, 
Lo) > eTA 4 x(0 dro), 





sc]. 5609) 


. ехр | lare s (eJ 9-9) + ars sei] 


+ ә(еј(—"-9)) J(-1-8) ) 


P(e 














-— 


- EXP E | are в(е)С"-9)) + arg p o9) 


To evaluate this, the following values should be used in this case: 


arg s (eX 6-9) A -A 
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arg s(sJ C9), “4 +2 
arg p(eJ 0-5) arg pled) 0 


pled) 22 


EI -— 


It follows then that 


XC — -6) 2 .2 ss, SE. 150 Is), e? 3 


~ 3 i (11-6 1( -Tr-Ö 
x e — 3 s (ef )) 222129 ә (ед à ) 
ər TT 
The desired value of the phase is about - 2” (or XE If the 
term in the brackets is greater than zero, this is what will result, 
If the term in the brackets is negative, the phase will be - = + TT 
=. ( or NOU T= -= он m. In other vords a branching 
2 2 Z 2 
will have occurred - not caused by the linear phase term of the 


groundwave, 


For the term in the brackets to be positive, it is required that 


s(eJ 67-9) > 9 gs(eJ T8) 


| A 
The value of 6 was such that 6 (n, +A, ) = 1/2 Oro > Ey e 
For a Fourier Transform of a discrete sequence, TT corresponds to 


N/2 in the DFT, This can be represented as a frequency Wo 50 that 
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W 


ô = - 
“0 
nun... O (in many applications it can be much larger), 6 ~ 20" 


What is required therefore is that 


"7 701 
s (eJ (53 н.Јӱ 0 (25-60 Wey 











Since |5(29") has even symmetry, what is required is that the 
Signal Fourier Transform magnitude decrease bysasraetor OL MOTE 
than 9 as w goes from .9875 Wa £071.0125 Woe The Loran-C signal, 
by itself, decreases by a factor of about 1.07 in this band. Consequent- 
ly, an analog filter which introduces an additional factor of about 
9 decrease is needed. 

lt can be seen that if the B4 used in the example were 0.9, 
the decrease would have had to be by a factor of — = 119), 


i5 n, vere greater than 40, the transition band would be narrower. 
Additionally, greater values of 6, i.e., 6 = ZA Eben 
will also cause this problem so that huge errors can accumulate. 

Before considering possible solutions to this problem, a few 
examples of the phenomenon should be presented, A 2048 pt FFT was 
computed for a simulated Loran-C pulse. The complex logarithm of 
the results are shown in the following figures, 

Figure B-13 shows an uncorrected phase curve for О 4 К < N/2 
with Ny = 0, no echo. Figure B-1/ shows the same for e 4, Figure 
B-15 shows the log magnitude and phase for ng 7 4, n = 20, B, = ,8, 
Notice that for values of k near N/2 (corresponding to w # m) there is 


a small amount of distortion. This represents the effects of aliasing 


as in figure B-10 - not enough to cause branchings. 
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Figure B-16 shows the log magnitude and phase curves for ng = ^, 
(n, uv» = 19,5. Here the problem can be seen (compare the last 
portion of the phase curve with figure B-6, Figure B-17 shows a corrected 
phase curve for Ny = 81: n, = 10, B, = 0.8, The phase unwrapping and 
linear phase term removal have been done properly. Figure B-18 shows 
the corrected phase for Ny = 0, (n,+4,) = 9,5, B4 = 0,8, The large 


error is apparent., 


A method of overcoming this problem is discussed in Appendix C, 
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Appendix C 
TRANSFORM PHASE CORRECTION METHOD 


This Appendix develops the technique used in this thesis to 
avoid the phase correction problems presented in Appendix B. 

One possible solution to the problems would be to use extra 
exponential weighting. As an example, if the B, in figure B-12 had 
a value of 1/3 instead of .8, the resulting vector representation 


Would be as shown in figure C-1, 


st ¿EN 
2/3 4/3 


A A A > 


Figure C-1 Vector Representation of P(e") With Extra 
Exponential Weighting 


In this case, the analog filter need only introduce a factor of 2 
decrease in the transition band. However, in view of the fact that 
the desired exponential weighting to be used would be a complicated 
function of signal parameters, and since the transition band would 


still be very narrow, an alternate solution might be desirable. 
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İn the cases considered in this thesis, the basic vaveform is 


knovn to be of the form pS 3 which has a phase value of almost 


- = for large w. From figure B-13 it can be verified that the 
phase atk x N/4 x 512 is essentially -1.57. In figure B-1 it can 


be seen that the phase at k = 512 is: 


arg X(N/4) s 1.7 - 44m x E 
y= = - 4(n/2) 


The scheme then becomes: examine the value of the phase at 
k z N/b, remove the -1.57 term which is due to the signal alone, and, 
determine the closest integer multiple of 1/2 in what remains. This 
integer is then ng. The phase can then be corrected by bringing 
the N/4 point of the phase up to zero by using the appropriate linear 
phase correction. This is done just as explained in Appendix B except 
that the straight line in figure B-2 goes from ihe k = 0 point to the 
k = N/4 point. The rest of the DFT record is then low pass filtered 
in a particular fashion - it is disregarded. A new N/2 point sequence 
is generated: 


0,1, ...,(N/4)-1 


|| 


Xə SƏ Ы Xa (Ə k 


xü (NES X-.C” DES Ee 


N 
new lat 2 Trees, (N/2)-1 


The effect of this operation is to double the effective sampling 


time (DT). A standard digital low pass filter is not used since zeroes 


jw 
in the sequence log X(e ) imply unity magnitude and zero phase in 
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the corresponding regions of the X(eI") sequence, This would produce 
an even larger error in the cepstrum. 

The effectiveness of this scheme can be seen by examining 
figure B-16 and noticing that the aliasing does not begin in the phase 
until well enough after k - N/l to be safe. A comparison of the 
resulting cepstra for the two methods of phase corrections is shown 


in figures C-2 and C-3, 





21524. 


3 Cepstrum for no echo: 
DT 2.2,0 - this becomes 
the à priori estimate 





4 | Cepstrum for one echo: 
DT = 2:0) Ta = 2020, 
0 
3 I, m 19 + 39754 B4 = 26 By 


minus a priori estimate 





Figure C-2 Cepstrum Computed With Conventional Phase Unwrapping 








DT = 5.0 - This becomes 
the a priori estimate 


EIUS 
Cepstrum for no echo: 










1 0 4 8 12 16 20 
Cepstrum for one echo: 
5 JO SIMON 3. 0:0 
Dİ - 5.0, B,İ- .8”B 
JA. ao 
minus a priori estimate 
-.5 


1 Cepstrum for one echo: 
» DT = 5.0, B, = +8 B, 


minus a priori estimate 





Cepstrum for one echo: 
To = 11.5, T4 = To + 40,0 


22 ОТ = 5,0, В, = .88 
| minus a priori estimate 





Figure C-3 Cepstrum Computed With Modified Phase Unwrapping 
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Appendix D 


MORE GENERAL NOISE DFT DERIVATION 


In this Appendix, the appropriateness of the white noise model 
of Chapter 5 is examined to see if the results depend critically 
upon it. 

Suppose the actual noise can be modelled as being obtained in 


the manner shown in figure D-1. 


H(t) Cao | ee) pi ei 


DT 


Figure D-1 Model of Non-white Noise Sample Generator 
Here, w(t) is a zero mean, complex, white gaussian noise process and 
h(t) is some appropriate function which can model the noise structure 


and any receiver filters. The analysis can proceed as follows. 
2 
E (+) vs) Ut 


E «| =0 


For the purposes of this analysis it will suffice to let 


A AO 
so that 


oo 


t 
f wy) h(t-y) dy = BS wly) EY) dy 


00 


v(t) 





537 


It follovs then that 


E vo) = 3 E İ «()] e* C dy € 0 
and 

E [v | SE | v(n-p2)| = 0 
Additionally, 


E с v*() | 


| 
TD 


1 f E [wQ) (y) PEN c - dy ал 
= g^ g* T f 6 (A-y) тб Оли-А-у) дуал 


For samples of the function, 


- nDT mDT 
R (nm) = E İv, vx|-8^ c [f J e) A a) 
00 =00 


For m> n this becomes 


n: DT 


R, (n,m) _ g^ 52 | ¿B(n-DT + meDT -2\) on 
32 2 ES DT + m DT - - TI 
2p -.00 


B g^ ,"8(n-n) DT 
2 


Similarly. Ior m< n, 
2 
R (n,m) = EL e78 (n-m): DT 
V 2 
so that 


к, (0) = E em. DI where 1 = |m - n | 
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The results of Chapter 5 depend heavily on an analysis of the 
variance of points in the DFT - from this, the noise magnitude bound 


is derived. Therefore it is important to consider R, ük) = E E jJ: 


a] | 


i -k -k (N-1) 
[ioe ма" oes "FM dəli 





O 1 2 


- 


z + vs v, + res Va cto... + и401-1) Һај 


-—: R (0) + REE m Tr ee ди X 


+ NO e? 3 i | 


+ es... 





Net Zerkm 
= N R (0) + 2 DI (N-n)-R (m) : COS N 
m=1 
N-1 
D-1 R (k.k) = Bo 15 x xr La" ə uu 
R m=1 


The evaluation of the sum in this expression can be simplified by 
making the approximation (N-m) x N. To see the validity of this 


approximation, the following would have to be specified, 


Bx 4.5 x 10° DT = 10 usec 


ane DT 


BD 7 and 7. 


The h(t) filter with such a pole location would have a half power point 





Fig. D-2 
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(N-m) 
“el “lu 
5.50 
A RA k=2 
k= 
kel 
-32 =16 I 
ly 8 TS 15. = 
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at a frequency of about 25 kHz so that p has a reasonable value, The 
sum in equation (D-1) then becomes: 

N-1 


EX (N-m) a" cos 
m=1 


Zııkm : 
N ş 





5. 


A plot of the three terms in the sum is shovn in figure D-2 for 
N = 256. Since the rapid decay of the a! term will dominate the sum, 
it can be seen that (N-m) Y N is a clearly justified simplification. 


The sum to be considered is then 





N-1 
N E a” cos em 
N 
=] 
: : 10 
which can be evaluated as in Jolley 
DE. — 1-.a cos 8 + P^ ББ GT BID ca co € 
2, a cosmo = > 
m=0 1 - 2a cos ô ft a 


IE E 296 and 299.2259 AL Az s 7 7. 50770 mis 


N-1 ə N-1 2 
20 ascos mi O = % a cos mö - i 
m=1 m=0 
2 
E a cos 0 - a E a [cos o - a] 
кк ә ERE - 
1 - 2a cos ð + Ld h -ae Jg 


The result is plotted in figure D-3 for 6 = ark 





o” 
24k Y a”) 


— Jemk|2 
N 





— 


1.45(80° 5 





Figure D-4 Non-White Noise DFT Variance 
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What can be noticed is that the zero crossing is at 


Q = oS 
or at N ud Zıık 
k = a; COS 223 = ¿214 N (since 9 = = ) 
It follows from equation (D-1) that 
2 N-1 
R, (k,k) ә Вә. N |1*2 02 a" cos m0 for a € .223 
^ m-1 ? әд: 


This is plotted in figure D-4. 

What can be seen is that as B:DT increases a decreases), the 
plot flattens out, approaching a value of N-R (0) as in the white 
noise case» The crossover point, i.e., the eae of k beyond which 


the noise variance is less than in the white noise case, was derived 


fron, 
OR. sf” 1 p 
Ko ı = cos E 
As B-DT—o, er DIET so that kọ approaches 57 2 = N/4, 


As described in Appendix C, it is the values of the DFT for k « N/A 
that are to be retained to compute the cepstrum, Additionally, the 
sampling rate for any given SNR is io be determined by ihe noise 
statistics in the DFT for values of k near N/!, What the analysis 
has shown therefore is that for kx N/4, the variance of the noise is 
nearly the same in both the white and non-white cases so that the model 


used in Chapter 5 is appropriate, 
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